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SUMMARY 

i  - 

,  *  :■> 


Homomorphic  deconvolution  is  applied  to  separating  the  com- 
ponents  of  a  convolution  in  which  multipath  Rayleigh-wave  propagation  is 
modeled  as  the  convolution  of  a  signal  S  with,  a  multipath,  operator  m.  The  * 
technique  is  applied  to  a  number  of  synthetic  waveforms  with  fairly  simple 
multisource  and/or  multipath  Rayleigh-wave  characteristics  to  demonstrate 

»  j  9  ' 

the  separation  .obtainable  between0S  and  m.  The  technique  is  also  applied  to  ; 

'  *  *  / 

several  events  recorded  at  LASA  and  ALPA  and’ results  of  the  method  as 

•  '  'V)  .  ' 

applied  to  actual  Rayleigh-wave  recordings  are  presented.  (  \ 
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RAYLEIGH  -WAVE  MULTIPATH  ANALYSIS 
.-USING  A  COMPLEX  CEPSTRfeM  TECHNIQUE 


SECTION"! 

INTRODUCTION 


The  excitation  and  prorogation  of  Rayleigh  waves  throughnilane 

»  '  "  -  i  '* 


multilayered' media  is  well  j^er stood.  Early  studies"of  Rayleigh-wave  record- 
ings  yielded  information  relating  to  the  average  crustal  structure  between  the  • 
epicenter  and  the  recording  site  by  group  velocity 'measurement.'  That  analysis 

.  *  **  .  *  5  ^  *  - 

allowed  definition  of  gross  structural  differences  -between  oceanic  and' conti¬ 
nental  crusts.  More  recently,  phase-velocity  measurements  have  been  Used  to 
determiner -crustal  properties  beneath  a  given  recording  sitei  If  the  "dispersed 

Rayleigh  wave  propagates  along  the  great  circle  path  from  epicenter  to  record- 

. 

-  ing  site, then  only  two  recording  stafions  at  that-  site  are  needed  for  phase 

.  *  'f  ^  ‘  *  .  -<  .1 _  •- 

velneity-determinaHim  '  The  -f  jcf-  that:  T?ayT<>Tgh  waves  do.  not  necessarily  pro- 

'  pagate  along-the  great  circle  path  but  may  b.e  laterally  refracted  through  highgr 

velocity  structure  and  refracted  at  continental  margins  which  ma'rk- the* boundary 
.  -  .  •-  ■"  »  *  -  _/ 
between  oceanic  and'continental  crusts  was  well  established  by  Evernden(19o3, 

-  '  *.  *  -  .  '  _  -  *  ' '  y  *. 

1954).  'A.tripartite  ar,ray  is  then  needed  to  a c cu raCely^determ i n e  the  direction 


■> 


-  form  techniques  for  calculating* phase* velocity  and  direction  of  approach  sub-. 
*<7  '  l  \  *  /  - 

stantiates  the  deviation  of  direction  of  approach  of  Rayleigh  -waves  from  the 

'  .  o  V  '  ♦  ■}  ~J  /'■  '  <* 

J  ,  *  ^  -  • 

great  circle  path  for  many  recorded  events  and  in  addition  shows  the-direction 

'j  T  •*  «  '  - ''" 

of  approach  to  be  frequency  dependent.  f  ,  _____ — - 


Lateral  refraction,  as  well' 'as  reflection  and  refraction  at  crustal 

*  4 

1  *  *  ‘  /  »  »  / 

1  transition  zones  along  different  propagation.paths,  yields/complex  interference 

„  _  i  ‘  *  -  _ ,  /  _  A 

patterns  of  the  recorded  Rayleigh  waves.  Recorded  v/aveforms^ofteri^apperr 


to  be  amplitude 'modulated  and  studies  by.  Pilant^and^Knopoff  (1904)  show  this 


may^ocgUT-because  ol-a  temporal  and/orspatial  distribution- at  tlie  source  of  / 


services 


\  A 


V  > 


7 


■  7 

f. 


Tf  I 


because  of  multipath/ 'propagation.  •  At  periods  which  corresponei^tobeats,  where  y  ^ 
.  a  minimum  o£  energy  is  observed  in  the  recording,  the  amplitude,  spectra  show  '  ,  | 

/  .*  *  . *  ,•  7  n  J- 

«  large  minima  and  irregularities  in.  the  smoothness  of  the'phase  spectra  occur/  tj 
These  phase  irregularities  produce  large  scatter  in  computed  phase  velocities  •'  .  ; 

.  — - *’ - ri— "  LI/  ; 


and  special  techniques  must^be  used  to  reduce  the  influence  of  beats* 

♦  Av  -  •  *  « 


• ..  i  - 


,  /f requency- wavenumber  analysis  by  Capon  .(1970)  yields  additional  f'/l  i 

*  ■  ■  •  (  ■  ■  y  j :  i  JJU 

evidence  of  multipath  propagation  which  is  frequency  dependent  and  conjectures  , 

were  made  concerning  the  actual  paths  taken  by  the  ^various  frequency  groups  j  r  /  : 

~  ~  ,  --  *  .  *  ^  ,  j  »  -  ,  •  l  ^  \ 

arriving , simultaneously  at  LASA  during  successive  200»sec  time  intervals'  over  e  '  ‘  ‘ 

'  v  '  5 

the  ‘Rayleigh -waveSijlu  ration.  ‘  '7  '  *  „  •'  *  V*t  ; 

*  "  -  . -  ■  '  ,  -  \  !•  ‘  i  •  ; 

/-  Measu'remehts'made  from  the  recorded  Rayleigh  waveform  to  i 

determine  group  velocity,  phase  velocity  and,  source  de  jth  are  strongly  influ-  .  •)  ;  j 

"  '  -  *  ,  .  /  ,  .  ^  .  < '  .  V’ 

enced  by  multipath  propagation  yielding  results  which  are  difficult  to  interpret.  In  ' 

•  < ‘  .  .  '  .  ■  -  •  ■  ’  I  -  '■«  .  " 

x  addition,  various  discriminants  calculated  from  the  contaminated  Raylei/gh  wav.e  |  ’  ' 

*  .  .  J)  WT  *  '  ...  ’ 

; _ (spectral  ratio.  'MS/AR,  RMS,  Mchirp,  etc. ),. are  not  representative  oy  the  / .  -  •  - 

//  •  '  .  _  .7.  -"o'"*"  -  ■  •  -  7  '<•  /  n  4 

Rayleigh' wave  excitation  at  the  source.  '  /  ^  •  /  /  |  t  | 


Rayleigh' wave  excitation  at  the 


source. 


Wl 


•  The  analysis  presented  here  then  is  concerned  with  recovering  /  /. 

'■  v  .  v  ii  I 

the  signal  and,  the  multipath  operator  from  the  recorded  Rayleigh  waveform.  /  « 

*  ,  ,  *  .  *  /  , 
Homomorphic  deconvolution  using  a  generalized  concept  of  linear  filtering  rs  }  :^ 

*  «i  1  ^  w  *  .  ^  “«/  1  t  1  * 

applied,  to  separating  the  components  of  a  convolution  in  which  multipath  pro-  l  *  ,  * 

,  pagation  is^modelbd^as' the  convolution  of  a  signal  with  a  multipath' operator.  ■  j-  7 


*L 


7  ,? 


The  signal  would  be  more,  representative  of  the  "pure"  Rayleigh  wave  and  the 
'^/kmiltipath operator' could  be  irfte rp reted  inj&rms-of  differential  dispersion 

J  -  -  J. _ -I- - —^=-7 

— and/ o rTqidtiple-sourcesT;  .  •  ■ 


■■  -  Lf 


i 


.  .  \ 
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SECTION  II 
TECHNIQUE 


,  '  Schafe'r’(1969)  has  presented  a^riew  approach  to  separating  convolved 

signals.  The  approach  attempts^separatiorTusing  a”complex  cepstrum."  techni- 

SJ  ■  -  ’  '  '/  /*  *  ‘ 

que*.  He  gives.'a  detailed  analysis  of  homomorphic  systems  for^decon  volution 
and  applies  the  technique  to.  a  class  of  signals  in  which  one  of  the  members  of 
the  convolution  is  an  impulse  train.  The*  technique' is  conceptually  simple  and 
Is  summqri'zed  as° follows:  ,  /  _  - _ „  _ - T-  "  \ 


The  problem  here  is  to  e stirffate^a^sighaTStt)  which  has  been 
complicated  by  a  few  .multipaths^m^tWthis  is  the'  multi- 
path  operator)!  .Mathematically. the  problem  is  X(t)  =  S(t) 
-*m(t),_  where  *  i^cb^vol^iti'on  and  X  (t)  is  the  observed 
waveform V/  ^ 


The^  or -Fourier  transform- is  taken  and  thus,  X(z)  =  S(z)  m(z) 

v’  :  :.//  ~  ’ 

TheTog~of^X(z). is  taken  4nd^the  imaginary  partms  made  conj 
tinuou s  and  thus^  log  X(z)  =  log  S(z)  +  log  m(z). 


Log  X(z)d.s^  treated  as  a  complex  time  series  and  a  linear  ■' 

-\  ... 

filter' is  used  to  separate^log  S(z)  and  log  m(z).  Bascially, 
thi|>re quires,  that  the  spectrum  of  log  S(,z)  and  log  m(z) 
ynot  ovdrl^ap- greatly.  «’  *  • 

The  exponential  is -taken  to  get  back  .to  an  estimate  of  t 
S(z)  or  m(z)'  There  is  no  ambiguity  here' as  exponentiation 
will  map  any  one  of  the  multivalued  log  function  values  b^ck 
irtto  the  same.cpmplex  nqrriber.  j 

Inverse  transform  ,to  get  an  estimate  of  Sft)  or  m(t).' 


.services  groyp 


This  technique  has  been  used  successfully  by"  Schafer  in  separating  speech* 
signals  into  their  more  b^sic  components  anctln-rieimoving  echoes  from' speech' 
data'.  The  cruciaMinkrn  this  technique  is  whether  or  not  log  S(z)  and  log  m(z) 
are  separable  for  multipath  Rayleigh-wave^afa.  -  The  numerical  procedu^e^ 
given  by  Schafer  were  programmed  (or  the  S/360  computer  and  synthetic  wave^ 

‘  J 

t  •  forms  representing  fairly  simple  multipat^i  Ravleigh^ave  Recordings  were 

used  to  determine  the  sepa>ration  obtainable.  The  sequence  of  numerical  calcu-,, 

•sx  .  '  ,  '  /'  •  '  -• 

lations  involv^tf  in. implementing  the  technique  is  described  in  the  following 


section. 


.  *  X 


The  sign  of  X(o)  is  checked  and  if  negative,the  si£»n  of  X(k)  is  changed  to  remove 
the  constant  phase  component.  Thi^-is  remembered  and  the  sign  of  Y(k)  is 


subsequently  changed  in  the  inverse  system.  The  complex  logarithm  is  taken 

'  ‘  \  '  '  . 


X(k)  ■=  log  [X(k)]  ?  log  jX(k)|'  +i  ARG  [X(k)]  '±  i2r  q  k=  0,1,...,  N-l 


and  the  principal  value  ARG  (X(k)']  when  q  =  0  is  obtained,  where 


-  ir  <ARG  [X(k)3  <r 


The  complex  logarithm  is,  multivalued  and  niust  be  defined  such' that  there  is 

‘  *  ,  '  ~ 

.  *  4  , 

no  ambiguity  with  respect  to  its  imaginary,  pact  by  choosing  the  integer  q  which 
makes  the  discrete  phase  values  approximate  a  continuous  phase  function.  This 
may  be  accomplished  even  for  rapidly  varying  phase  angles  provided  the  fre¬ 
quency  sampling  is  fine  enough  by  either  one  of  the  two  subroutines  SCHAFR  or 
AARON'r  SCHAFR- uses,  a  first  difference  and  AARON  uses  a  one-point  prediction 

"  .  l 

error  to  detect  jumps  in  phase- of  2*-  .  The  in^put  signal  is  augmented  with  zeros 

„  .  *  v  *  > 

to  give  an  input  trace  N  sample  points  in  length  in  order  to'  obtain  the  fine  fre¬ 
quency  sampling  needed.  The  value  N  is  a  variable  input  parameter  to  the  pro-, 

*  * 

gram.  The  estimated  phase  curve  is  then  rotated  by  removing  the  linear  phase 
component  and  arg  [Y]is  subsequently  reverse  rotated  in  the  inverse  system. 

"i  *■ 

Next  the  inverse  fast  Fourier  transform  (IFFT)  is  taken. 


N-l 


i  Zw  kn  / 


x(n) 


l  y*  X(k)  e  N 

v'  \ 


N 


k=0 


./ 


n  =  0,  1, ; . . ,  N-l 


\ 


-A  - - 


-yielding  the-  co_mplex“ce~pstrum  x. 

'  ‘  *  V  . 

The  parb  due  to  S(t)  in  the  complex  cepstrum  is  primarily  concentrat- 

.  ,  ,  •  i 

ed  in  the  "short  time"  region  while  the  impulses  attributed  to  m(t)  occupy1  , 


the  "long  time"  region.  If  these  signalsjare^  separated  in  "time",  then  a  linear 


filter  system  can  be  used  to  separate  the  signals  from  one  another.  A  short- 
pass  or  comb  filter  system  can  be  used  to  remove  the  impulses  due  to  m(t) 


from  the  complex  cepstrum  while  retaining  most  of  the  part  attributed  to  S(t). 


i 


Similarly,  a  long-pass  filter  system  can  be  used  to  remove  most  S(t)  while 


i 


.  The  filter, systems  are  applied  in  the  program  by  zeroing  the  appro¬ 
priate  portions  of  the  complex  cepstrum.  The  sharp  cutoff  of  these  filter  systems 

•  ’  A 

will  produce,  ripples  in  the  transform  of  y therefore  an  option  is  provided  to 

modify  the  filter  amplitude  function  with  a  cosine  taper  giving  smoother  transi-^,  • 

1 

tions  between  one  and  zero.  .  ,  ,5 

j 

o 

The  filtered  complex  cejjstrum  y  then  is  input  to  the  inverse  system 

,\  *  '  '  x  /  .  ,  .  - - —  — 

where  the  FFT  is  taken,  .the  imaginary  part  reyerse  rotated,  the  complex 

exponential  taken,and  the  IFFT  taken  to  yield  the  filtered  output  y.  < 

\  _ ^ 

Preprocessing  of  the  input  signal  can  also  be  performed  by  the 

•  |  *  ;  i 

program  when  necessary..  First,  the  mean  is  removed  and  the  signal'is  trun- 
cated  at  the  nearest  zero  crossing  at  each  end.  A  zero-phase  bandpass  filter 
is  designed  and  convolved  with  the  signal  to  allow  decimation,  the  cutoff  fre- 

/  .  .  •  /  i 

quencies  and  decimation  interval  being  variable  input  parameters.  If  the  i  __ — it 


input  sequence  is  nonminimum-phase, then  exponential  weighting  may  be  applied  :'v 
in  order  to  make  the  sequence  more  like  a  minimum  phase  sequence.  It  has 
been  empirically  observed  that  this  will  insure  Jhat  the  part  of  the  cojmplex  cep-  ’ 
strupn  attributed  to  m  will  occur  in  the  region  t  ^  t^  where  t^  is  the  spacing 
between  the  fir^t  and  second  impulses. .  The  program  applies  a  truncated  ex¬ 
ponential  window  to  the^input  sequence.)  ’  \  / 


;  ;  t  - 

w(t)  =  a  }c(t) 


0  s  t  <  L 

o 

elsewhere 


where  j  a  J  <  1  »  L  is  the  length  of  the  input  windo^r  and  a  is  a  variable 


input  parameter.  Subsequent  calculations  are  performed  on  the  exponentially 

weighted  input  sequence, and  the  output  of  the  inverse  system  is  then  unweighted 
N  '  ■  i.  ‘  ■  '  ' 

using  the  same  value  of  a.:  .  ,  /  \ 


In' addition,  the  program  generates  CalComp  plots  of  various 


functions  obtained  during  the  sequence  of  calculations  to  allow  visual  analysis  of 

I  ■  \ 

the  results.  "  • 
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has  bees  made  continuous  and*  then  rotated  to  remove  the  linear  component. 

The  IFFT  of  the  reaTand  imaginary  parts  gives  the  complex.cepstruim  For 
-  '  r  " 
this  example  the  input  X(t-)  is  minimum- phase,  after  rotation,  yielding  a 

•>  '  *  __  •  "5  , 

complex  cepstrum  that ds  zero  for  t.<  0.-  Th^  part  due  to  S^(t)  is  primarily 

•in  the  "short,  time"  re'gion  from  t  =  0  to  t  =  80.  The  impulses. due  to  the  .. 

■  '  &  •  ,  *  \  ’  .  , 

echo. are. in  the  "long  time"  region  and  occur  at  t  =  80,  t  =  160,  etc.  The 

>  complex  cepstrum  is  then  short-pass  filtered  at  ^sample  point  74  and  the  FF1I:~^'^ 

.taken'.  '  The  real  part  shows  the  effect  of  removing  the  rapidly  varying  component 
■  ‘  •  -  -  ■  -  ..  .  ..  ,  .  • ' 
due  to  m(t).  The  imaginary  parf^is  reverse' rotated,  the  complex^exponential  .* 

■  ■  .  ■  /  ••  \  *  .  .  .  -  J  -  •  • . 

taken  and  the  IF-FT  ..taken  to  yield  the  short-pass  gutput  which  is  the  estimate 

*  <  *.  ,  . . ,  •  .  ■  ,  ->  *■ 

of  S  (t).  An  error  trace  is  calculated  fS  (t)  minus  the  short>pass  output]  which 
o,-  'shows  th,at  the  output  of  the  short-pctSs,  system  estimates  S^(t)  j^ith  negligible 

.error.  Next,  a  long-pass  filter  at  sample  point  74. is  applied  to  the  complex  ,  '  , 

*  "  ]■*  •  .  4  o  -  .  » 

cepstfum  and  processed  as.  bdforeT  The  real  part  shows  thex  effect  pf  removing 

y  **  /  *  '  V  " 

the.  slowly  varying  component  due  to  S  (t).  The  output^of  th(b  lOng-pasS  system 

/  '  j  •  . 

-  '(bottom  trace,  ^Figure  IV- 1)  yields  mj(trwith  negligible  error.  •  .  : 


1 
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r  .  SECTION  V  %  -  -S. 

CEPSTRUM  ANALYSIS  OF  SYNTHETIC  MULTIPATH 

RAYLEIGH  WAVES  ’  "**  - 


Seven  syntactic  waveforms  with  simple  multisource  and/or  multi - 
?  -  __  --  -  . 
path  Rayleigh  wave  characteristics  were  processed  using  th&  complex-cepstrum 

/technique  in  order  to  determine  the  separation  obtainable  between  S(t)  and  m(t) 

-  ' 

and  to  serve  as  an  aid  in  interpreting  subsequent  outputs  obtained  from  actual 

Rayleigh  wave  recordings.  .  •  « 

*  a  _ 

9>  .  /  ■  X 

.  '  The  waveforms  used'were-produced  by  adding  together  two  or  more 

cosine-'taperfed  chirp  wav-eVorms.  ^  First”  a  chirp  waveform  is  generated  by 

'  -  ' ,  r  -  '  ■v 

specifying  the  initial  frequency^i^/'final  frequency  f ana  time  duration  L  , 
where  the  impulse' response  is  '  --  . 

*  .  r  -  /  -  %  .  • 

''*"/•  r  *  » ■  -  f  -  f  _  v  •  " 


.  ./  -  r  /  - 

P(tJ  =  sin  .2iT  [f  ■* 


r  “  '  *0 

V/  * 


)•] 


O^t  L 


‘  •’S.  otherwise 


where 


%  -fg  =  0.025  Hz  and  f  ^  -  C  ~”P  Hz  . 


#  *  O'  / 

Shown  in  Figure  V-l  (top)  -is  a  chirp  waveform  (R(t)with  a  time  duration  of 


500  sec. 


:■  -  ^  ^ 

In  order  to  make  the  waveform  more  nearly  resemble  an  actual 

■-  '  c  * ' .  * 

seismogranS  recording  of  a  Rayleigh  wave,  a  cosine  taper  is  applied  to  R(t)  where 

7  *  '  f  „  -  I  *  "  ,  . 

the- output  is  *  ‘f.  <  "  ' 

.  .  S(t)  =  R(t)  (l  -  <=<>=  ~TT")j  '  :  ,  0StS  L  • 


This  is  analogous  to  shaping  the  spectrum  with  a  recording  system’s  respons^/ 

’  v.  * 

(long-perio.d  seismometer)/  The  cosine-tapered  chirp  waveform  (Figure  V-l, 

*  -  ‘  '  ■  /  ‘ 
middle) 'resembles  an  actual  teieseismic  seismogram  recording  of  a  Rayleigh 

wave  in  that  it  exhibits  normal  dispersion  and  its '-amplitude  spectrum  is  highly 
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CHIRP  WAVEFORM 

At 3 1  SEC 

.  '  #* 

s** 

r~'  _  ; 


— tAPEREDGHiRP 
At  =■  1  SEC  , 

\ 


DECIMATED 


X 


^  — - : - - - - - — - : - ’ - — - 

•peaked;  '  '  .  •  •  . 

*  •  -  •  ** 
Application- of  the  cepstrum  technique  was  found  empirically ,  •. 

-  o 

yto  work  better  by  resampling  the  data  at  8-sec  intervals  (Figure  V-l,  *  bottom). 

This  sample  period  gives  a  Nyquist  frequency  of  0. 0625  Hz  w'hich  is.  well  out-. 

\  '  * —  /  .  \ 

side" the  frequeifcy  band  specified -for  the  chirp  waveforms.  As.noted  in  Figure 
V-l, resampling  of  the  waveform  yields  .an  irregular  appearing  time  function, 

\  *  -  '  J  *  *  ^ 

however,  there  is  no  loss  of  spectral  information  in  the  band  of  interest.  a  -Para--  ^ 
meters  calculated  directly  from  the  decimated’ time  function  such  as  group  -  . 

• 1  velocity  (wave  frequency  and  "arrival  times  of  pfeaks  a‘nd  troughs)  and  maximum  / 
amplitudes  of  given  peaks  and' troughs  are  determined  using  a  three-point  quad-  / 
ratic  interpolation  yielding  results  essentially  identical  with  those  calculated  / 

^  from  the  waveform  sampled  at  i  sec,.  ,  .  .  \ 

The  ffrst  multichirp  waveform  analyzed  consists  M  a  500  sec  chirp 
with  onset  at  1-sec  plus  a  500-sec  chirp  with  onset  at  250  sec.  Shown  in  Figure  V- 
are  the  two  individual  chirps,  their  sum, .and  the  decimated  sum  which  is  the 
_1  input  signal  to  the  complex- cepstrum  a:  lalysik  program.  '  Group  velocities  cal-  , 
cuiated  from  the  waveforms  are  presented  in  the.  right  side  of  the  figure.  For  . 

refeffence.the  group  velocity;  curve  (sol  id  line) -for  the  LASA  perturbed  model  _ 

’  :  ■  I .  .  i  ' 

(Texas  Instruments  Incorporated,  1.9674  is  included  in  the  plots.  A  distance  ^ 

\  '  }  '  1  *  t  " 

"  and  origin  time  is  chosen ‘such  that  the  resultant  group  velocity  pointy  calculated 

from  S  (t)\closely  fit  fhe  LASA  perturbed  curve.  For  thip  example,  the  distance 
*  i  ‘  V  \  *  * 

was  9053  Ion  and  the  origin  time  was^2543  sec  before  trace' onset.  7  The  multi - 

J  ,  •  '  y  ,  J  ■  '  \ 

chirp  waveform  yields  a  com.plicp.ted  group  velocity  curve.  ^The  smooth^  points  at 
at  low  frequencies  and  the  smooth  points  ait  high  frequencies  derive  from  the  non- 

\  *•  - * —  *  •  :  \  ‘  \  y  //  I  */  (- 

o\^er- lapping  portions  of  S^t)  and  S2(t)  respectively ./  A  complex  intexference 
pattern  results,where  S^t)  and  S2<t)  overlap  in  time/ yielding  the  scattering  of 


LASA  are  presented  in  Figure  V-3.  LASA  recordings  of  Rayleigh  waves  from 
events  occurring  i<n  mapy  widely  separated  regions  yield  complicated  group 
velocity  curves  with  characteristics  very  similar  to  the  multichirp  examples 


points  observed  at' intermediate  frequencies..  For  comparison,'  group  velocity 
•‘calculations  for  several  earthquakes  (catalogued  in  Table  VI- 1)  recorded  at 


presented  ip  this!  section. 
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.Figure  V-3. 

Group  Velocity  vs 

Frequency  Plots 
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The  multichirp  waveforms  are  not  minimum  phase,therefore  ex¬ 
ponential  weighting  is  applied  to  the  input  signal  and  subsequent  calculations 

f  , 

are  performed  on^the  exponentially  weighted  sequence.  The  value  for  a  used 

in  the  weighting  was  0.  98  for  this  example,  and  the  output  of  the  system  is  then 

* 

unweighted  usirig  the  same  value  foi  a. .  • 

*  \  .  '  .  , 

Results  of  the  sequence  of  calculations  using  the  complex- cepstrum 

technique  are  presented  in  the  remaining  portion  of  Figure  V-2.  Shown  are  the„ 

!  .  ■  .  . 

real  and  imaginary  parts  of  the  complex  logarithm*  of  the,  FFT  of  the  exponen¬ 


tially  weighted  input  signal.  The  slowly  varying  component  due  to  S.and  the  rapid 

\.  a  •*  •  j  j 

ly  varying  component  due  to  m  are  clearly  seen  injthe  plot  of  the  real  part  al* 

I,  ,  *  *  &  '  ’  ^ 

the  log.  The  value  N  given  in  the  figures  is  the  number  of  samples  for  the  FFT 
ancf  is  a  power  of  two.  The  input  signal  is;  augmented  with  zeros-before  trans¬ 
formation  to  give  a  tjrace  N-saynples  .long.  For  (the  frequency-domain  plots,N 
is  the  number  of  frequency  samples  and  the  folding  frequnecy  w  is  a  N/2.  'Fo;r 
a  sample  period  At  =  8  sec,  the  folding  frequency  is  •  ‘  ' 


=  '0.  0625  Hz  . 


and  for  the  time  domain  plots  (complex  cepstrum) 


N/2  Ns  NAt/2  =  819?  sec 


The  complex  cepstrum  is  shown  next.  The  impulses  attributed  to  m  are  small 

in  amplitude  relative  to  the  pulse  attributed  to  S  which  occurs  primarily  in  the 

\  ,,  \  -  1 

"sh^rt  time?’  re\gion;therefore  an  enlarged  plot  of  the  complex  cepstrum  is 

also  shown.  Tbe  large  pulse  about  zero  time  is  clipped  injthe  enlarged  plot 


for  display  purposes.  The  impulses  attributed  to  m  occur  at' t  =  249, 
t  ’=  498,  etc.  and  are  easily  discerned  in  the  enlarged  plot. 


~~  \  The  complex  cepstrum  was  short-pass  filtered  at  208  sec,  the 

'  \  \  0  I 

second  TFT  taken,  the  imaginary  part  reverse  rotated,  the  complex  exponential 
\  ■  ’ 
taken,  the  second  I  FFT  t^aken,  and  the  sequence  unweighted  to  yield  the  short-pass 
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Win, 
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i  f; 


'  output  which  is  the  estimate  of  S  (t).  The  error  trace  shows  that  the  output 

\  ‘  \  .  '  ■  1  ' 

of  the  short-pass  system  estimates  S^(t)  with  little  error.  Group  velocity  qalcu- 

*  «*»  *  i  ’ 

lated  frpm  the  short-pass  output  is  shown  to  the  right  of  the  time  trace.  - 


Next4:he  cor 


was  lone-pa  si 


processed  as  before.  The- Teal  part  of  the  second  FFT  shows  the  effect  of  re- 

,  '  l,  -  / 

moving  the  slowly  varying  component  due  to  S.  The  output  of  the  long-pass, 

<  ,  '  j  i 

system  yields  m(t)  which  is.  essentially  j  ,{t)  +  $( t-249)  for  this  example. 

■  ^  The  second  multichirp  waveform/ analyzed,  consists  of  a  600-sec 

♦  v  ‘  A  /  \  i  .  i  / 

chirp  with  onset  at  1  secTplus  a  600-s'ec  chirp  with  onset  at  250  sec.  Results 

/  '  .  — - -w - '  ■  . 

of  the  sequence  of  calculations  are  shown  in  Figure  V-4.  The  format  of  the 


plots  in  this  figure  and  succeeding  figures  in  this  section  follows  »that  of  the  , 

1  ’  ; '  \  ’  •  ■  ! 

•  first  example  given.  S  (t)  and  S_(t)  are\each  100  sec  longer  and  overlap  100 

r  •  i/ L  ’  \  ' 

sec-more  in  time  than  for  the  first  example.  The  short-pass  output  yields  S^j(t) 


with  practically  no  error  and.as  seen  in  the  figure, the  long-pass  output^rp(t) 
is  essentially  identical  L/ith  the  first  exaigp^e.  \  -  , 


x  ,  The  third  multichirp  waveform  analyzed  consists  of  a  500-sec-chi-rp- 

wi^h  onset  at  1  sec  plus  L  600-sec  chirp  with  onset  at  250  sec.  The  analysis 

results  are  presented  inIFigure  V-5.  The  short-pass  output  again  yields  S^(t) 

with  negligible  error.  The  long-pa^s  output  ni(t)  is  more  complex  than  the  first 

two  examples  where  a  constant  time  delay  existed  between  S  (ty  and  S  (t)  for  all 
•  \  1  I  ^ 

-  frequencies.  For  this  example,m(t)  can  be  interpreted  in  terms  lof  differential  \ 

V  ,  '  '  4  •  -I  \ 

dispersion  where,in  the  time  delay  between  S  (t)  and  S  (t)  is  249  sec~kt  0.  025  Hz 
\  „  ,  \  ■  1  .  *■ 
and  increased  with  increasing  frequency  to  349  s.ec  at  0.  050  Hz. 

The  fourth  multichirp  waveform  analyzed  is  more  complex  in  that 
1  /  •  \  _ 
it  consists  of  three  chirp  waveforms  ,500  sec, x 550  sec  and  600  sec  in  length 

\  •  V 

with  onsets  at  1  sec,  201  sed\and401  sec.  The  analysis  results  are  presented 

in  Figure  V-6.  The  complex  \cep£trum  was  short-pass  filtered  at  160  sec  and 

~~the  output  yields  S ^ (t)  with/ittle  error.  The  long-pass  output  m(t)  is  more 

complex  than  the  preceding  example  in  that  it  shows  the  differential  dispersion 

between  S,(t)  and  S_(t)  as  well  as  betweer^\S  (t)  and  S  (t).  As  note£  in  the  plot 
1  2  (  \  13  ' 
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Figure  V-4.  Cepst^urri  Results  for  Synthetic  Multipath  Rayleigh  Waves.  Input 

Is  the  Sum  of  Two  Chirp^Waveforms;  600  Se,c  and  600 *Sec  in  liength 
with  Onsets  at  1  Sec  and  250  Sec  and  Amplitudes  of  1  Unit  and  T Unit 
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\Figure  V-6.  Cepstrum  Results  ;for  Synthetic  Multipath  Rayleigh  Waves.  ‘Input 

t<  4  *  Is  the  Sum  of  Three  Chirp  Waveforms,;  500  Sec,  '550,  Sec  and  6 00  Sec 
'  in  Length  with  Onsets  zg  1  Seo,  201’  Spc  and  401  SevC  and  Amplitudes 
of’l  Unit.  1  Unit  £nd  1  Unit.  !■'  ‘  h'  i' 


the  time  delay  between  S,(t)  and  .is  200  sec  at '0.025  ,H z>  and  increases  wit 

-  ,  '  i '  2 .  i *  x-  ;  j  .  •  ' 

increasing  frequency  tc^ 250  sec  at  0.  05.0  Hzwhile  the  time  delayybetw^een  S ^  (t) 


and  S^tyis  -400  s.ec  at  0.  025  Hz  and  increases  wjth  rising  frequency  to  500\^  * 

sec  *at  0;  050  Hz.'- •  -  v  '  \  ’ s>  ^  ^S'X>.  .  \ 

•  ;  \ - '  '  '  • 

-  •  •  ■*  .  ■  '  i  \  * '  “K‘\  c.  . 

,  ;  ‘  'The  fifth  mulfjchirp  waveform  analyzed  is>identical  to  th€~first 

■  -  -  *  *  ‘  \  *  -  /V  * .  f 

example  with  the.  following  exception.  The  waveform's.  (t)‘ is  scaled  to  one- 

-•»  .  ■  '  t‘  ,  I  ..  •  ;  •  .  -  '  • 

half  the- amplitude  ofjS^t).  This. is  clearly  reflecte’d.Th  the  lon^-^jass  .putput 

for  this  example  (Figure  V-7)  where  mft)  is  essentially  J  (t)  ~’\6.5$  (t-249). 

\.  -  -  - T,  •  '*  •  .  '*  \ 

'/  /  .  \  - 

i  ’  -_  The‘ sixth  multichirp  waveform  analyzed  Consists\qf  two  chirp  wave-  • 

—  "  11  *  r'  i.mi— -  —  —  .  |  r—  tm  ”  *  %  •  ^ 

,  /  *  ••  '  '  V  *  ■*  »  '  f  4  ’  '  »  * 

forms  500  sec  in- 'length  with  onsets- at  1  sec  and  101  'sec'  ( Figure.  V- 8,)..  .The  time 

»\  .  ’  f  '** 

separation  .between  St  (t)  and.S  (t)  is  much  less, than  in.  previous. example's,  tn 
/  - ,  i  • "  •  1  .■*>  2  .  j-  *  «  "  .v  - 

addition,  S  (t)  is,  scaled  to'on'e-half.the  amplitude  of  S'  (fc)  and  the  Waveform  in-  - 
2  '  '  '  6  -  •  1  ’  J  .  « 

'1^  '  |  ,  ,m  .  '  m  »  ‘  •  “V  , 

.verted' with  respect  to  Sj  (th  The.  complex  cepstrum  is  long-pass  filtered  at 

'  *  *  *  *  *  .  .  A  •  ^  % 

80  sec,yieldingih(t)  which  is.  essentially  5 -(t)  -0. 5  5  (t-100j,*  -The  copnplex 

j  ^  •  •  ’  _  !  i 

cepstrum  is  also  . short-pass  filtere&^'at  80  ;s‘e  (^yielding  the  output0 £ h own The-  - 
calculated  error, for  this  output;  although,  still,  fairly  Small,  -is  largpr-than  ob-  •  . 

'*  „  '  |  .  *  o  s  -  '  •  »  „  *  • 

,  •  .  '  ^  ^  'V  •  *  ,  ■  -» 

tained  from  the  previous  example's.  The  snort-pass  system  has  removed  all 

^  ,  *  #  -  *  *  ,  ,  _ _ 

of  the  contribution  causedby  m  and  part  of  that  caused  by~S.  :  ;  *  f- 

-  ’  -  \  i  ,  .  f  i'4  '  '  ’  • 

\  [v  *  *,  ,  *  ^  ,  »** 

'Nl\"  ,  •  y  ^  '  %  4  °  ,  , 

.  The  sey^nth  multichirp  waveform  analyzed  c opsists-aLteO-chirp  wave 
forms  400  sec  and  600  sec' in  length  with  onsets  at  1  sec  'and-h'sec The  analysis 

*•'■’[  .  \  ;  <  x  *  t  ‘  ' 

results  a,re  presented  in  Figure,  V-9.  For:.this.. example  the, time  delay  between 

'  .  o  !  i  ’  -T  |  *  ^  ■  •  '  -  t  ^  i  , 

S.  (i)  and  S_(tJ  is  zero  ht  O.J325'Hz  and  increases  with  increasing  frequency  to_ 

^  1  "  2  ■  f  .  '.'-■s',  ,  , 

200  sec-at  0.  050  Jlz.  An  {inspection  of  the  cOmplex  cepstrum  shows  that  S  and 

-/U-  ..' .  «  .  .  ..■  X' 


s^t)  and 
200  sec-i 


m~ overlap  si 

.attributed  to 

*  ' '  '  ■  V* 

to  S  and'wice 


pb^'^ec  waS  applied  to  .the  complex  cepst.rum.  The  filter  at  56  sec  yielded  the 

*“'*>!*  *  *  *  \  /*  '  ■ 

smallest  error.trace.  -  The.shQrt- and  fbng-pass  filtered  results'  for  that  filter 

■-  f  .  -  '  *,  j  '  n  1  O  - 

are  shown  in  the  figure.’  ,  Simple  short-  or  lon^?1pass  filtering"cahnot  completely 

’  •  i  ^  »  '  ?  ,  / 

separate^  from  rp.  Comb  filterii-.g  could  be  Used  to  remove  the  part  attributed  _ 

1  ^  1*4  7  *  '  *  /  -  /V 
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Cepstrum  Results  for  Synthetic  Multipath  Rayleigh  Waves*  Input 
Is  the  Sum  of  Two  -Chirp  Waveforms;  500  Sec  and  500  Sec  in  Length 

with  Onsets  at  1  Sec  and'250  Sec  an4  Amplitudes  of  1  Unit' and  1/2  Unit 
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V-R.  Cep. strum  Results  for  Synthetic  Multipath  Rayleigh  V/ aves.  Input 

Is  the  Sum  of  Two  Chirp  Waveforms;  500  Sec  and  500  Sec  in  Length 
with  Onsets  at  1  Sec  and  101  Sec  and  Amplitudes  of  1.  Unit  and  1/2 
Unit.  (The  Second  Chirp  Wavdform,  with  Onset  at  101  Sec,  is 
Inverted  with  Respect  to  the  First  Chirp  Waveform. ) 
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Fioure  V-9a  Cepstrum-  Results  for  Synthetic  Multipath  Rayleigh  Waves.  Input 
°* .  \  /  Is  the  Sum  of  Two  Chirp  Waveforms;  400  Sec  and  600  Sec 'in  Length 

'  with  Onsets  at  1  Sec  and  1  Sec  and  Amplitudes  of  1  Unit  and  1  Unit 
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to  m  while  retaining  almost  all  of  the  part  from  S,provided  m  was  a  fairly 

J  *  '  ■  '  1  ’  \  / 

simple  impulse  train.  However,'  for  this  example,m  represents  a  transfer 

l 

*  v  ,  *  *  - 

function  of  considerable  length  with  zero  time  delay  at  the  ‘lowest  frequency 

(O' 025  Hz).  The  part  attributed  to  m  in  the  complex  ceprtrum  is  a  complicated 

» 

series  of  overlapping  transfer  functions  occupying  the  "short  tinge"  region 

*  to  which  comb  filtering  cannot  be  successfully  applied  When  the  two  components 

i '  *  i  *  ° 

significantly  overlap,:tKen  the  ..linear  filter  system  removes  too  much  of  the 

_ _ -r-  '  '  *  4  *  ;  '  •.  .  _ . _ 

complex  cepstrum  of  the  de’sired  output.  However,  for  this  example,  the  signal 

.  .  ,  •  /  * 

estimate  obtained  is  clearly  a  better  approximation  of  S^(t)  than  is  the  contamina 

ted  input  signal  and  discriminants  (AR,  RMS,  chirp  statistic,  etc. )  calculated 

I  ■  I  I 

9 ^  from  tpe  signal  estiinate  will  more  nearly  approximate  the  true  values.  In  • 

addition,  it  should  be  pointed  out  that  the  group  velocity  points  obtained  from  th£ 

/  *  ■  (  •*  ' . i. 


short-pass  output  (Figure  Y-9)  closely  fit  the  expected  group  velocity  curve  \ 
■'  of  S  (t).  /  •  •  •/ 

.  The  results  presented  in  this  section  demonstrate  the  feasibility  of 

this  technique  for  fairly^;simple  multisource  and/or  multipath  waveforms.  Thi 
multichirp  inputs  tested  represent  simulations  of  teleseismically  recorded;  ' 
Rayleigh  waves  resulting  from:  ,  *  1 

-  Interference  from  two  or  more  seismic  sources 

occurring  at  the  same  location  with  small  .time 
separation  and  no  multipath  propagation  .  '  ' 

!  <~V  i  \ 

•  *  Interference  from  sources  at  different  locations 

which  arrive  at  the  recording  site  with  small;  •  ~ 

A  /  ~  “time 'separation  *  * 

•  Interference  from  a  single  source  generated  ' .  \ 

Rayleigh  wave  which  propagates  along  twa  or  \  \ 

•  three  distinct  paths  with  differing  group  velocity  \  ^ 

f  characteristics  therefore  arriving  at  the  recording  \ 

site  with  small  time  .separations 


The  examples  show  that  for  cases  where  $  and  m  are  separated  i^. 
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"time',’  a  linear  filter  system  can  be  used  effectively  to  separate  the  two  com- 

x  '  _ 

ponents  from  one  another.  The  technique  requires  no  a  priori  estimate , of  Sland. 


\ 

will  give  an>  estimate  of  both  the  signal  S  and  the  multipath  operator  m. 
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SECTION  VI 


PRELIMINARY  CEPSTRUM  ANALYSIS  OF  RAYLEIGH 
WAVES  RECORDED  AT  LASA.  AND  ALPA 


Several -events  recorded  at  LASA  and  ALP^were  processed  using 

the  complex  cepstrum  technique  in  order  to  make  a  preliminary  evaluation  of 

_  «■ 

the  method  as  applied  to  actual  Rayleigh  wave  recordings.  The  events  used  had 
good  signal-to-noise  ratio  and  associated  information  for  the  events  (location, 
date,  origin  time,  etc. )  is  given  in  Table  VI- 1.  The  vertical  component  of  the 
Rayleigh  wave  was  time  partitioned  from  each  event  recording  for  input  to  the 
program.  >  \  \ 

. . . •  . '  '  \"  r  ’ 

Processing  results  for  event  NTS03  recorded  at  ALPA  are  presented 
in  Figure  VI- Shown  in  the  figure  are  the  decimated  input  signal,  thei  real  and 
imaginary  parts  of  the  log,  the.  complex-cepstrum,  the  shoi^-pass  filtered  out¬ 
put,  and  the  long-pass  filtered  output.  The  input  signal  has  a  modulated. character 
and  minima  are  observed  in  the  plot  of  the  real  part  of  the  log.  These  Vminima 

could  of  course  result  from  either  an  interference  phenomenon  or  low  signal 
.  ’  \ 
excitation  at  the  source  at  those  frequencies.  Whether  or  not  energy  at  those 

frequencies  was  present  at  the  source  cannot  be  ascertained  from  this  beafned 

output  alone,  however  an  inspection  of  the  group  velocity  curve  calculated  fi*^m 

this  event  (Figure  VT-2,  top)  shows  that  the  recorded  waveform  is  not  a  simple 

normally  dispersed  Rayleigh  wave  and  this  characteristic  is  indicative  of  inter- ( 

fering  waveforms. _  •  , 

0  x  * 

The  complex  cepstrum  fof  this  event  is  much  more  complicated  than 

-T*  \ 

those  obtained  for  the^simple  examples  given  in  the  previous  section.  The  part 

of  the  complex  cepstrum  attributed  to  m  is  not  a  simple  impulse  train.  Thus  it 

would" be  very  difficult  to  apply  a  comb  filter  system  since  this  requires  much 

information  about  the  number,  locations  and  durations  of  the  impulses  in  the 

complex  cepstrum.  An  alternative  is  to  use  a  short-pass  system.  A  suite_of 

# 

'  .  '  "  ‘  1 

short-and  long-pass  filters  covering  the  interval  from  -  16  sec  to  —  128  sec 


29 
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was  applied  to  the  complex  cepstrum  and  outputs  obtained.  An  analysis, of  those 
results  was  made  by  inspection  of  the  plots  ^showing  the  real  part  of  thJ  sWcond 


'  FFT  after  short-and  long-pass  filtering.  -The  short-pass  filters  at  16  and  20 
sec  removed  all  of  the  rapidly  varying  component,  due  to  m  and  .also  removed 

•  ’I  |  \  . 1  *  , 

much  of  the  part  due  to  S.  The  filters  at  36!  to  1128  sec  passed  almfost  all  of 
the  component  due  to  S  and  also  passed  increasing  amounts’  of  the  .rapidly  ' 

•  .  I  \  ■  .  .  *  ’  , 

varying  component  due  to  m  for  increasing  filtdr  times.  The  filters  at  24‘ to  32> 

•  c\'  v  ,  * u  \  '  [  A  \ .  ■  .  .  •  ,  .  < 

sec  passed  most  of  the  part  .due  to  S  ana  removed  almost  all  of  the  Rapidly  vary¬ 
ing  component  /  Shown  in  the  figure  are  the  reW  part  of  the  second-  FFT  and  the 
output  for.  both  short-and  long-pass  filters^af^JS  sec. 

,  ■  *  •'  v  •  A  •  *  *  *4  ■ 

vy  \  *  ^ 

*  .  •  i  t  \  ; 

I .  Group  velocity  was  calculated  from  the  short-pass  output!  A  compari-, 

/  •  ,  \  <  ,  v  i  '  , 

son  between  the  group  velocity  points  obtained  from  the  input  waveform  ■and'the*  : 

Ishort-pass  output  (Figure  VI-2,  top  and  bottom  respectively)  shows  that  the* 

ol  /  *  \  '  \  r  \. 

iscattef  of  points  is  largely  reduced  for  the  output  trace  and  .the.  points  are  much  .  , 

■  v\  -\  ;  j  ■  1 

rore  continuous.  'In  addition,  the  points  closely  fit  the  theoretical  curve  calcu- 

/  '  *  *  "  *\  - im¬ 

itated1  for  a  continental  crust.  The  great,  circle  path  for  this  event-  traverses  , 

l  j ,  ,  j  *  i  -  •  • 

Nevada,  Oregon,  Washington,  parallels  the  Coast, Mountains  of  British  Columbia, 
a^id  crosses  a  portion  of  Alaska.-  '  ,  '  ,  !  ^ ^ 


1/  The  iOpg-pasS' output  show?  a  large  complex  impulse  at  about  30  sec. 

interpretation  of  the’  long-pass  output  based  on  several  assumptions  could 
oftcpurse  be  made:  However,  additional  information; would  be  requij^d  to  sub- 

'  if  v'-'  .  •  *  -  *  .  i 

stamtiate  that  the  assumed  conditions  and  therefore  the  interpretatior^were 

,  correct.  J\or  example,  two  sources  with  identical  excitation  and  l’CTcation 


collect.  Fpr  example,  two  sources  with  identical  excitation  and  l’CTcation  « 

occurring  with  small-time  separation  w&uld  yield  a  long-pass  output  as  Shown 

1  ‘  /  tr  \  "  '  / 

*  \  •,  %  U  1  \ 

fo,r  examples  one’and  two  of  the  previous  section.  However,  it  is  conceivable 

|  -  '  y  .  }  ‘  .  ■<  J  '  •> 

thatlan  entirely  different  physical  situation  Oould  yield  an  essentially  identical,* 

wavlform  recording  and  subsequent  long-pass  output.  Therefore  an  analysTsr'\^  . 

of  snort-period  p  waves  recorded  at  a  network  of  stations  and  subsequent  hypos 

center  calculations  would  be  necessary  in  order  to  confirm  or  deny  the  possibi- 

lity  that  two  sources  of  the  same  type  and  magnitude  occurred  at  the  same  location 
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.  ^  j  /  #  - .  *  '  &  *  ■ 

with"  the  time  separation  observed,  in  the,  long-pas's 'output.  Conversely,  mow-  - 

•\  -  v‘  I  ;%  i 

'  ledge  of  the  information  described  above  would  allow  accurate  filter  placement 
in  the\complex-cepstrum  technique  tOiyield  thL  optimum  signal’  estimate x 

,  - "  '.\  '■  -  ’  .  f  ■/  ;  • .  -  ;v  • 

*  „  .  '  <.Next,  a  Greenland  Se:a  event  L-01  recorded  at.LASA  \tfas  professed. 

•  \  .  .-..•>!»!  1  ••  ; 

Those  results  are  presented  in  Figure  VI -dJ.  The  great  circle  path  for  mis  event 

*  \  '**/';*  \  f  jjjg  ■  * 

traverses. a  portion  of 'the  Greenland  Sed,  ’Greenland  and  the  Canadian  SHijpld. 

Group  velocities  calculated  from  the  input  signal ’and  the  shprt-pass  out™*  a^e 

\  —  ■  •  .  |  r-  1  jjf  <  /.  4, 

[ shown  in  Figure  VI-4.  The  input  signal  is  relatively  free  of  beats  and  is&a.  fair-' 

!  v  ■  -  \  •  .•  a  -  .  4  '.r  /  :  -  ■  J  fife, 

1«r  o  i  nrtvmolltr  ^  C  "fiOPn  in  flip  Cf  Tflllll  Vpl  fin  fv  W  '  H  GW — 


ly^simple  normally  dispersed\waveform  as  seen  in  the  groug  velocity .  ^ How 

>  ’  .  -  ‘  ,  ,  #  \’  ’  .  *  "  ,  iMr..  i  ‘ 

ever/  interference  is  occurring  late  in  the  waveform  and  a  large  scatters  pre-r 
t  .  \  *  •  <-  -  ■  •  *■  l  iff/  .  h  \ 


'  'n.  r  -  >  \  .  *  •  .*  I  ‘  '  •  *  *)  &  ■  \  \  .  ‘ 

sent  at  the  low  group  velocities  which  corresportd  to  the  coda  of, the-  event.  •  Group 

.1  .  ‘  »•*»  .*•*  '  V*  It « .  ,r  *  \  * 

\  ^  0  t  »  *1  J  .  I  -  f  * 

velocity  for  the  short-pass  output  is  somewhat!  smoother, and  as  seen  iMFigure  VI- 3, 

’  \.  *  ’■  !  I  j  1  1  ■  .\ 

much  of  the  coda  of  tjie  input  event  is.  not.  pre*  vrit  in  the  signal  estimate*  The  long- 
pass  output  shows  a  large  complex  impulse  -at  about  60  sec]  The  longj ringing 
train  around  -300  to  600,  .sec  is  attributed  tip  the  most  rapidly  varying  .component, 

*  ■  •  *  j  '  .  *  •  '  ;  j 

i.  seed  id  the  plot  . of  the  real  part  of  the  log,  with  a  Af  of.  about  0.  002  tp  0.  003  Hz. 

5  .  '  «  .  /  *  * ,  “  I  •  '  *11!'  \- 

»  0  j  “  ^  ^  ‘  1  /  -  * 

Four  events,  from ’the  Kurile'Ts./  recorded  at  LASA  weref  processed.  . 

*  >■  .  \  ’  ».  \  ■,  ' ,  j  /  .  •  X-,-  Ti  V 

Rayleigh  wave's  record.ediaLLASA_ff.om  events  . in  this  region  are  riotpriously.com- 

I  /*.'T  /  ■ 

plex,  and  parameter s-calculated'  frbm  the,  waveforms  give  results  which  are  diffi¬ 
cult  to  interpret  as  evidenced  by  the  group  velocity  plots  shown  in  Section  V„  Fig- 

’  I  s'  /  •  ‘  , 

/  ure  V-3,  The  great- circle' path  traverses  the  Kurile  Arc,  Bering  Sea,  Gulf  of  .. 

;•  '  i  J  ^ 

-  Alaska^  and  Western  Canada  —  a  corifplex  series  of  structural  environments. 

. .  r  ■ 

Complete  cepstrum  processing  results  for  LL-2008  and  LL-2018^  are  shown  in 

i  '»•./  ■  •  *  _ _ L/ 

Figures  VI- 5  through  VI- 8.  .  For  brevity, ''’only  the  input  signals,  “shortrand  long- 

pass  outputs,  and  calculated  group  velocities  are  shown  for  events  LL-2019  and 

- - - ; _ .  ^  —  -I 

LL-2028  (Figure  VI- 9).  All  four  events  show  beats  in  the  input  signal  and  .  ,  _ 

—  .  *■ 

^minima  are  observed  in  the  plots  showing  the  real  part  .of  the  log.  Again  the 
complex  cep^tra  are  not  simple  impulse  trains  and  are-cSfficult  to  interpret. 

A  suite  of  short-pass  filters  was  applied  to  the  complex  cepstra  and  outputs'  S 

calculated  as  before.  Presented^  the  figures  are  the  outputs  for  a  filter  at 
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igiire  IV-4.-"  Grojup  Velocity  v§  Frequency  for  L-01.  (Top,  Input 
'  Signal?  Bottom,  Short-Pass  Output) 
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Figure  VIS.  Cepstrum  Results  for  Event  EL- £008  Recorded  at.LASA 
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igure  VI-8. #  Group  Velocity  Vs  Frequency  for  LL-Z018.  (Top,  Input 
■  _  ’  Signal;  Bottom,  Short-Pass  Output) 


80  sec.  i  he  signal  estimates  obtained  are  clearly  much  simplier.  than  the  input 
•  * 

waveforms.  The  LL  event  recordings  .obtained  from  Lincoln  Laboratories  were 
bandpass  filtered  at  0.025-0.055  Hz  during  preprocessing.  Since  there  is  no 
energy  above  0.  0625  Hz,tne  event  recordings  were  decimated  by  8  giving  a 

V 

sample  period  of  8  'sec  'and  a  Xyauist  frequency  of  0.  <3625  Kz  at  input  to  the 

complex  cepstmxn  program.  Group  velocity  points  calculated  from  the  input, 

waveform  sampled  at  8  sec  showed  variations  with  respect  to  the  waveform  . 

sampled  at  1  sec.  Those  variations  occurred  only  at  the  high'er  frequencies 
.  '  '  ,  *  * 
where  frequency  is  changing  very  slowly  with  time  for  these  events.  The  three- 

point  quadratic  interpolation  used  in  the  group  velocity  calculations  to  determine  . 

*  » 

peak  md  trough  times  is  adequate  over  the  frequency*b'and  of  interest  for  wave¬ 


forms  sampled  as  coarse  as  4  sec  and  no  smoothing  was  performed  on  those 


1 


waveforms  wi thSt  £  4  sec,’  howeve.r  for  the  waveforms  sampled  at  8  sec 

V.  m 

smoothing  of  the  order  number  vs  time  function  over  about  5  points  was  necessary 
to  reduce  the  variation  in  calculated  group  velocity  at  the  higher  frequencies.  It 
should  be  pointed  out  that  smoothing  over  5  points  does  not  significantly  alter 

the  group  velocity. points  calculated  from  the  original  waveform  sampled  at  1  sec. 

* 

'  '  •  n  * 

Group  velocities  calculated  from  the  input  signals  and  the  short-  , 

pass  filtered  outputs  are  presented  in  Figures  VI-6,  VI-8  and  VI-9.  The  simpli- 
*  *  *  » 

citv  of  the  group  velocity  curves  for  the  outputs  is  apparent..  The  curves  obtain- 
.  -  -  *  •  **.*-*.,* 

ed  are  neither  oceanic  nor  continental  but  are  an  intermediate  variant  of  the 

*  C  '  * 

two.  Group  velocity  calculations  are  based  on  the  epicentral  distance  between 

>  ' 

event  and  recording  site  which  assumes  a  great  circle  path  for  the  surface 
waves.  If  the  first  arriving  surface  waves  (at  a  given  frequency)  have  devia.ed 
from  the  great  circle  path  by  traversing  a  higher  velocity  structure, then  the  cal¬ 


culated  group  velocities  will  be  lower  than  for  the  actual  path  taken. 


The  events  were  long-pass- filtered  at  80  sec,yielding  a  fairly, 
sm.ple  output  for  several  of  the  events  wherein  4  or  5  separate  complex  impulse 


trains  may  be  observed,  ihe  similarity  of  the  outputs  obtained  is  strongly 
suggestive. of  several  fairly  distinct  propagation  paths  from  source  to  receiver  , 


A? 


sery toes  group 


/ 


for  these  events.  This  type  of  analysis  coupled  with  frequency -wavenumber 
analysis  (Capon,  1970)  could  provide  a  powerful  tool  for  understanding  the  com¬ 
plex  nature  of  recorded  Rayleigh  waves. 

•  -  * 

*  * 

Six  additional  event  recordings  were  processed  using  the  complex- 
cepstrum  technique.  The  results  of  these  calculations  are  presented  simply 
by  showing  the  group  velocities  calculated  from  the  input  signals  and  those 
calculated  from  the  output  signals  (Figures  VI-10,  VI-11).  The  results  obtained 
are  similar  to  previous  results  in  that  the  outputs  are  simple  normally  disper* 
wavetrains  as  evidenced  in  the  group  velocity  plots.  _ v 
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Figure  VI-1'1.  Group  Velocities  for  Input  Signal  and  Short-Pass  Outputs. 

Filter  at  64  sec.  Input  Signals  are  Single-Channel  Record¬ 
ings  of  Event  EPX  14646 
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APPENDIX  A 


COMPLEX  CEPSTRUM  PROGRAM  DOCUMENTATION 


I*  DESCRIPTION  *  1  ^ — ■ 

Homomorphic  deconvolution  using  a  generalized  concept  of  linear 

filtering  is  applied  to  separating  the  components  of  a  convolution.  The  method 

'  . 

is  based  on  transforming  a  convolution  of  waveforms  into  a  sum,  using  a  linear 

’ 

filter  system  to  separate  the  .additive  components  and  transforming  the  result 

...  .  '  '  1  .  " -  .  . 

back  to  the  original  input  space.  An  extensive  presentation  -of  the  approach  is 


, given  by  Ronald  W.  Schafer 


(19<?9),  1 


Echo  Removal  by  Discrete  Generalized  - 


linear  Filtering:  Technical  Report  466,  Massachusetts  Institute  of  Technology, 

‘  *»  .  -  k  _  .  1  I  V 


Cambridge, 'Mass. ,  jPhD'Thesis.  .  * 


tyP19a^  1:1111  time  for  the  complex  cepstrum  program  required 


approximately.  3  minutes. onjhe  360/65. 

H  .  •  -  -  *  - 


obtained  during  the.  sequene 
sis  of  the  results. 


-  -  *  "  '  ,  ’  t  * 

CalComp  plots  bf  the  inputs  and  outputs^  as  well  as^various Timctior 

/  *  •>  / 

j  the.  sequence  of  calculations  are  generated  to  .allow  visual  anaiy-  • 

?  »  **  1 
c.  ?  I 


II  RESTRICTIONS  .  *  . 

ft  .  '  * 

.1)  The  input  data  trace  may  not  exceed  1946  points.  ,  , 

2)  The  program  is  dimensioned  for  a  maximum  of  2048  complex' 
—  ^ 

■  nurhbers  in  the  FFT  calcelatipr.s. 


HI  INPUT  «. 
Card  1- 
•  Card  2 
Column 
1-10 


Card  3 


Card  4 

1-10 
11-1 5 
16-20 


NN100 


Variable 
NKEPS>  ' 


NAME 


IDIN 

NSKIP 

NPTS 


l  ;* 

»  *  *  * 

Mode  Description^  ’•  v  ^  ' 

/ 

1  •  Number  of  event,s- tg  process  in 

’  this  run.  //  *  ■ 


A,‘  '  »  AloHanum^ric  title 


;'\nput  Trace  I>E^  • 
P\  ’-<j[  .^Number  of  point 


.^Number  of  points  to  skip  on  input 
>  Number  of  points  to  read  on  .input  • 


services  group 


s 


<M  (M 


V*  *■  -  *  ^ 

Card  4  (continued) 


■f  * 

•V  -  V  ' 


* 

,  Column 

.Variable 

Mode 

Description  .  -  jy 

r 

21-25  ■ 

* I XL  -  - 

if 

I  * 

%.  **  9 

Decimation  interval-  after  band  limiting 

26-30 

!  FI'LT 

>  • 

'  I  _ 

=0  stop-afte"r  calculating  the  complex*  T 

• 

» 

cepslrym.  -  .  - 

*  '  ' 

. 

=>0  number  of  filters  to  apply  - 

31-35  ... 

-I  AARON 

I 

■  =1  call  Aaron’s  algorithm  **  ’ 

*  „  °  *,5  ,R 

Jr  1 

=2  call  Schafer’s  algorithm 

£ 

9 

*  • 

(for  correcting  phase)  '  .= 

36-40 

-  XEEX  '* 

I  ' 

Yota!  data  .Je?hgtf»  for  transforms,  etc. 

0  m 

• 

(must  be.a  power  of  two).  Zeros 

m  *  •  « 

*  « 

are  added  to  give  an  input  trace  '  % 

m 

- 

NLEN  points  in  length. 

..'41  -,50  -■ 

DBLIM 

* 

r. 

.  >0  clip  real  part  of  log  at.  DBLIM 

db  below'peak.  « 

*  «?. 

<0  no  clipping.  '  „ 

51-60 

ALPHA 

E 

=  ‘  a.  Exponentially  weight  -the  -input 

• 

« 

trace  X.  'W(ri)  =  a n  X(n).  .  • 

•6.1-70 

*  DT 

*  E 

,  Sample  period  of  the  input  data.  \ 

"  c 

Card  5 

*  * 

■  '  '  *  •  \1 

lrlO 

IDOUT . 

,1  ■- 

Qutput  ID  for  data  to_be  saved  on  tajie^ 

11-20 

IDIS 

,  I 

Output  disposition  for  data  .to  be-  saved 

- 

*  v 

on-tape.  ' 

/Card  6 

’  %  _ 

Filter  car.ds,  Supply  IFILT 

*  j 

times  *  ■’ 

1-TO 

ITYPE 

T 

▲ 

=  1  short-pass  system. 

•  » 

• 

=2  long-pass  system.. 

t 

v  * 

« 

=  3  comb  filter  .  ♦ 

1 1  -2a- 

I  NIT  PT  # 

.1 

Point  at  which  to  filter  the  complex* 

• " 

cepstrum  (if  ITYPE  =1  or  2).- 
Number  of  comb  filter  points 

N 

' 

(if  ITYPE=3).  t 

2.1-30 

IDX  ' 

I  * ' 

2*lDX-fl  points  will  be  zeroed  at  each 

*  c 

•  - 

-  •  .  comb  location  (supply  if  ITY‘PE=3). 

‘31-40 

^ — '  JTAPER 

I  ' 

•  ^O'cosine  taper  the  long-or  short- 

V 

< 

pass  filter.  .  1  » 

$* 

=  1  no  taper  on  filter. 

Card  7 

• 

Comb  Filter  cards, 

(Format  8110)  supply  if  ITYPE=3). 

\  Supply 

(INOTCH(J),  J^l,  INITPf)  comb  centers. 

3 


A-2 
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rto<Ziot\*  „ 

T  i-/>cr  r^>!  or.<>rP<>oT<>fPFrOFILTE9<>FFT<>ROT  QEXPOIFFTOJT 
THIS  Jc,  4  TFGHMIOIIF  DEVFLOPMFNT  VERSION  WHICH  CAN  BE  SIMPLIFIED. 
APDP.priABLY  AS  MFl-'JL.  AS  SPEEDED-  Up  FOR  PRODUCTION)  TYPE  RUNS 
.  0!'*FNSin\,s0  Fno  ?043  REAL  "DATA  POINTS  *  •  *  . 

DIVISION*  ARFAIBOOOl.NAMFI 15).AUM(7|  J  *  ' 

.DIVISION  *  I  20^8  I « WORK  C  20491  ,  %  ' 

DIvENSIPN  rPMPI?04R|  .  '  « 

OIMFNSfPN  PORI  2048) 

DIMENSION  S"AVr  1-204  9 )  .  '  »  : 

DIMENSION  T5ACcjf?049l  •  *  .  .  . 

HI  ,S.T0RI?p4R*  '  *  •  . 

.OI«f'NSIr’N  I  NOTCH ( 1.00 J  ’  .  '  v 

o(Mpf.'«;n»j  TAPFpIBI  -  ,  _  V' 

,  ,-s.ypLFX,  X.  SAV=  *  -  . 

DATA  f  APFR_y.  0948.  . «397,  .9660 ».  7660*  .6424-*. 5000 ,.3420*..  1736/  c 

*\f  ALL  NNI00  -  .  ’ 

1  FDOM  AT  ( 1£.  I°RE  1  S.-7 1  '  ,  ’ 

p  I  =’.141593  ?  \  "  > 

.  TWPPI  =  4. 093194  ’  -  , 

‘CALL  .PI  OTS  I  AREA. 9000.3. 01  •  #  -  '*  *' 

CALL  DAY! IOaVFJ  ‘  '  '  ‘ 

»EA0.C,.7)  NKFP5  '  ■  '  '*  .  -  ‘  .  t 

Dr  10nn/i<FPS=l  .NKFPS  #  .  .  ' 

PALI  STASK  •  •  '  . 

'  RFADI5.  I  )  'NA«  .  " 

1  FOOMAT  f  1 5A4I  '• 

WP  I  T  F  (  ft.,?  )  NAME  . 

7  FO°M  ATl  I  HI  ,  1  X  ,**NAME  ’  =  *.15A4)  -  ^  . 

W°  I T  c'f.fi .  3 )  I  DATf  •  '«  *  ' 

3  pro  M/1T.I  /./  »  l  Y  .  *  OATF  =  •  *  1 1 0  J 

PF  AD  15’.  4  ) ,  I-PIN.  NSK  I P,  NPTS,  TjNC.IFILT,  I  AARON.NLEN.  DBLI M*  ALPHA,  OT 

4  F  OP  MAT  (  I  19.6I5.3F10.0_J_  f  ‘  ■  ■ 

9  FO°mAT  (  /  /,  IX.  •  I  DIN.  NSK  I  P..NPT5,  I NC  ,  I  F  ILt  ,  l  AARON, NLFN.  D8L  I M,  ALPHA, 

W°  I T  c (  5. 5  )  I  0  IN.N-SK  I p,  NpTS»  INC  »  I'F  ILT,  I  AARON,  NLPN-.  DBLI-W,  ALDHA  ,DT 
9  P0°M  AT  {  1  X,  3[  10.  F.10-.2.-E  I0.4/1F10.  2) 

of  AD  (5, 7)  IDOUT.iniS.  .  '  '  ,  . 

7  FOR  “AT  f  >  I  1  0  1  ’  ‘  !; 

WRITFI6.9)  inOUT.IDIS  r  ...  '  .  *  ' 

q  prpM/\r  I//.1X.'  I.DfJU  T  ~  t  II 0.'  •,  10  IS  =  *»I5) 

0FLY-7D.49/f-Nl  FN*INC*OTJ  *  „  ” 

■  -  f \<f  =0cc IMAT I CN  R  AT  F  •  - 

DT  =  SAMP|  F  pFPIpO  OF  INPUT  DATA 
'  NPTS=NUMRFR,' pF  I NP.UT ''PQI NTS  ’  ^  ’ 

NPF=10J  ’  '■ 

FOrP.n  ’ 

C.  l>r  .  O'*  r<  '  ‘  - - '"'S  •  *  '  , 


services  group 


.V 


r 


k 

t. 


is  «*i' 


E2-.oso  "  -  * 

F3=.0 f>?5:  .  . 

WRI-TFI6.13)  F0.E1.F2.E3 

FORMAT  (//*  I Xf  *  CUTOFF  FREQUENCI FS  FOR'  PREF ILTFR.*  *  4F10.4) ' 
DESIGN  *7 FRO  «HASP  BAND  PASS  ‘FILTER  :  ' 


t //« IX**  i.UH.JI-h  f-H 

DESIGN  *7 FRO  fHASP  BAND 
FO^l  j/  f  FLOAT*  KPF.9*DT). 
FP  =  l./(?.*OTl) 


470 
‘  475. 
‘490 


-  NF=l+INTfJF9/5FQI  •  ; 

FTFST=F0*F4  ^  '  **..'»* 

F=0.  -  •' 

no‘  480  I_=1,NF,  -  *'  *'•  /  '  .  ■. 

‘IF  rFTF.ST  1460**460*440  -  *  -  •  -  ' 

'  I F ( FO-F )  490*449*445'  *  *r“T'  * 

HOPK  CD  =0*  .  '  .  *  .  -  ' 

GO  TO  *480'*  .  J 

i  IF(FItF).  460.455*455  .  1  y  ./•  v  ' 

;  wp«.Km  =  n.+cns.(Pi*cF-En/«o-E'i/m/2,.  I 

•Gfi -TO  490'.  •>  * 

I  IFCF2-F)*  470,465.465V,  ?  /  *  . 

WIRK (!)»!•  »  .  •  V  l  ^ 

GO  TO  480  •  c 

IFIF3-FI  445,475.475  -  *  ' 

.  H0RK(T.l  =  U.-C9SCPI*CFr63)/(.E2-E3m/2.  -  ^ 

i  F=F*FD  '  •  * 

LF=(NPF-1 )/2  ..  *  •  ' 

.V=-OT*FinAT(LFI.  '  *~  '  ?  .  . 

no  490  I =1 .  NPF 

TEMPm=HbRK(l)  .  ‘ 

F?FD  •  '  *  - 

bn  485  J  =  2*NF ' 

TEMP!  I  F=TEMP(  'H+2.^C0S(T«0PI.*F*YI*H0RK(dl  1  * 

I  F-'F+FD  X/  *  ‘ 

YiV+OT  ,  / 

I  TEMPI T l=TFMPf I |*FD*OT  .  ‘ 

READ  TR  ACE »  TRUNCATE  FACH  END. .  AT 0 NE AR EST  Z  ERO  CROSSING 
APP.I.Y  COSTNF  TAPER  TO  EACH .  END','  BAND  PASS  FILTER.- 
DECIMATE.  RFMOVE  mean 
CALL  RDCPI IDT N.WpRK.NSKIP.NPTS.l.AUM) 

AMAX=0. 

SUM=0.  '  ' 

DO  ’406  T  =  1 ,NPTS  * 

.  SUM=SUM+WnRK ( I ) 

SOM=  <  SUM /NpTS ) 

nn  407'  I  =  1  ,  N P T.9  '■ 

T  WORK (I  I =HORK f I i-SOM  '  *•  ' 

INDEX=l  „  •  '  -  t. 

on  815  J=1.NPTS  ;  ' 

IF (HOPK (J) .GE.O. • AND. WORK! J+1J.LE.0.)  GO  TO  816  .  ‘  • 

I F  (  WOR  KlJ).LF.O.  .  AND.  WORK  (  J+l)*GE.O.)  GO  TO  81*6 
5  I NDF  X=  I NOE  X+-1 

•  GO- TO  1 009  v  . 

i  TF  f ABSCWORKI INDEX)) .GE. A8S( WORK ( I NDEX+ l ) ) ) INDEX? INDEX+1 
.NPTS=NRTS-INn£X-fl  .  „  A 

no  «17  J  =  1 » N  P.T S  -  ,  "  . 


sarvicas  group 


i 


I 


\ 


«i  7  work  (j)=*orkUnoex-i*ji. 

*  INDF,X=.NS<  f  £>♦  INDEX- 1  t  . 

WD  I.TF  {  6  *  49,1  )  INDEX  '  *  ' 

4Q I  FORMAT  (//.  1*.  •  TOTAL  NUMBER  OF, POINTS  SKIPPED  ON  INPUT  =**.1101 
tN0F^=*iPT<  '  V  ‘  . 

J=1  , 'I P  TS-  -V  •  .  ‘  f  - 

„  fFJ  wn£>K(INOc<)  .OF.O..AND.WORKC  INOEX-ll.LF.O.)  HO  TO  819% 

I F  { tfno  v  {  index)  „LF .Q.  .'AND.  WOR  K I  I NDEX-1 )  •  GE.9-.1 )  GO  TO  819  ' 

818  1  NDF  x:==J  NOP.X- 1  .  '  ‘  "  ... 

GO .  TO  1^8  *  '  ’ 

819  ir^l  ARS  (WORK!  INOcX))  .GF.ARS1  WORM  I NDE  X-rl  I  >  >  ,1  Nf)EX=  INDEX-1  . 

*JPT S- 1 NDF  X  S  -  •  '  . 

|-*Jt>Tc  -7f  '*  *  - 

;  '  ‘  K  =  8  * 

on  870  -1  .8 

'  WP°K  C.fl  -  f  APFR  (  11  ♦WOPRl  J  »• 

/  w'PR  K  ( K  1  -  T  A  PF  8 ,1  1  1  *  WnRK  C.K  > 

’  .l  =  /J  ♦  1  .  '  % 

8?0  K=<-  ]' 

NHA  F=NPF  h?  •  • 

00  4?8  *1=  1  «NHAF 
478  TRACE!  I  )=0.. 

'DO  8  78  <J=1  ,NP-TS  ’  ■" 

878  .TRACE!  NHAF»'j)=WGRK(J  )  ,  ■"/ 

on  ioo  j=l,NPTS  * 

mo  WP^K  ( J 1  =  TRACE ! NHAF+J )  -  ’ 

CAL L  WRCf  P  !-I  mtlT  ,  WORK  ,NPTS .  1 . 1 , I DIS .  IDATE,  NAME,  AUH) 

„  •  .  ind)UT=TDOllTM  ’  '  -  ,  ' 

k=nhafojpjs>,i  *  • 

.  l.  =  K  *NHAF-t  .  ‘  *  ’  . 

np  4  80  T=K.I 

480  TRACF!  I  1=0.  .  ■* 

r  <=PI  l  T.PR  INPUT  TRACE 

?FRn=TFMPf NHAP^l 1 

no  4,81  1=1  .*NPT<;  '  •  ,  . 

“  <  4  *■»  , 

"  .1=1-1  \ 

‘  K  =  NPFf.l  ...  '  .  .  ‘  . 

'SUM  =  0..  '  •  l  ^ 

■  D(1  478  J=1  .NHAF  '  .  '  * 

476  SUM=Sl»M^TFMP(  J  )*(  TRACE  (  J  +  Ll+TR  ACE  ( K-J))  '  ' 

481  WORK  (1  )  =  SUm*7FR0.*TRACE  INHAFM)  1 , 

r  *  D  F  C  f M  A  T,  F  THF  BAND  I.  I  HI  TFT)  TRACE’ 

1=0  .  ,  ”  •  '  \  . 

DO  408  1  =  1  .NptS,  INC  '  .  .  ,v 

>  1=1*1 

408  TP  AC  F 1 l  1=WC8K ( 1)  .  4 

NP  =  l  . 

/■  PUNCH  THF  RANT  L  I  M  j  TFO»  DECIMATED  INPUT  TRACE »  FORMAT  (  5E  1 8. 7  ) 

CALL  WPCP(  I  opUT  ,  TR  AC  F  tNP  *  l.,  •!  ,80,  I DATF  »  NAME  »  AIJM ) 

I  P  T  S  =  Np:*  1  *  .  *  ‘  '  ’ 

LPTS=NP 

NPT  S=Nl EN  , 

r  -  p x ° AND  TIME  SEALE  ON  INP^jJT  DATA  AND  OUTPUT  DATA 


/ 


FX  PA  NO  =  4  . 


A-6 


services  group 


4 


KPTS=\»PTS/F*PANO  “  ( 

AEtX=DFL x'***PAm  *'  •  -  .  '  .  ‘  : 

AT=DT*  INC  v‘  '•  '*  . 

WR I  TFI6,V01 1  AT  ...  *  .  ■  -  „ 

40  1  FORMAT*//:,  IX,*  SAMPLE  PERIOD  OF.  BAND  LIMITED  AND  DECIMATED  DATA* 
IF  10^,21  f  '  ,  *  ,  ~  •*  .  • 

ho  700  J- T PT S'* NPTS  ♦,  ■  ?.  ’  -  .  "\ 

700  TRACEl J )=AMAX  ’  .  *  '  .  '  '  .  ' 

r  '/  PI  OT  INPOT  TRA#E  '  ’  ’  ;  t. 

CAL!  LETTER^ b.»0.»0.2*NAME*90.,60)  '  *  .  '  ' 

YSIZF=2.  -  .  -•  "  ‘  ‘  U  -  ^  .  • 

IN=0  •  " '  .  *  *  ;  - 

!30X=?  -  .  '*.■ 

XORIG^YSTZF+1 1  <  .  '  '  .  ■ 

CALL  ORIGJNIX0RTC.O.I  ‘  ’  -s.  ■ 

CALI  LE7TEPI0.75.0., 0.15, ‘INPUT  TRACE*  ,90.,,11)  * 

CALL  PLOTA I  TP  ACE»KPf  S»'AELX,  YSI ZE.*  XLEN*YMAX  ,YMI N,iN, I  BOX  INC,  DT ) 
CALL  WPXPC1D0UT, TRACE, NPTS,1, l.!DIS,IOATE, NAME, AUM) 
!DOUT=IOOilT*l  •  '  .  J  •  '  .  V*  -  , 

■  thax=yMax  **  •  #  f 

.  TMIN=YMIN  . 

DO  7*50  J=l.LPTS  , 

.  Aj-j  — .  •  fy  t  *  ’ 

750  MQRK  ( J!=TRACei JI*(ALPHA»*AJI  / 

00  7-70  J=fPTS,NPTS  .  x  ; 

770  WORK! J»=AMAX  «  ■ 

C  .  c  STORE* REAL  INPUT  INTO  COMPLEX  ARRAY  X  ' 

DO  8  J-l ,NPTS •  -  / 

8  XI J|=CMPLX(WORIC  l J) ,0.01  ^ 

C  TAKE.  FFT  OF  ARRAY  X 

CALL  NLOGN I X, NPTS* — 1 • 0 I  ’  . 

C  TAKF  l  OG.  OF  X  '  '  ,  . 

-TFNO=NPTS/?+l  v.  *  •  , 

ISW=1  .  <•'  ’  /  ■ 

- RF AL  P=RFAL( X 1 1  1 1 

TFIRFAIP.LT.O. I  I  SW=2  * 

GO  TO  CllO, 1201, TSW  ,  . 

'  170  no  no  J=1,IEND  • 

RFALP=-REAL(XC  JM  '  ‘  '•  ♦ 

,  AT=-AIMAG!Xt J)t  * 
no  XI  JI=CMPLX(RFA10,A!)  • 

WRITFIA.1A0)  .  . 

l'AO  FORMAT  i  IX,  ’CONSTANT  PHASE  COMPONENT  REMOVED*) 

1:10  CONTINUE 

00  11  J=1 , 1  END  -  * 

11  XI 4) =C LOG  I  X I J) )  ‘  •  ' 

RFAL.P=RFAL  I X I  I  END)  1  ‘ 

.  X 1 1 FND. )  =CMPl  X{i?/:ALP,0.  )x 

rmax=reac'(xi  d  )  <•  •• 

TNnEXP=l.  ‘ 

DO  17  j=l  v-TEND 

.  REALP=R6AL!XIJ)1  *  .  .  > 

WapX{.l)=PFALP  j  3  *  ,  . 

ifi.realp.lt .rmax  )  Wm  1? 


services  group 

s  ■ 


- 


RWAV=«PAL-P 

INOFXP=J. 

];?  continue 

'pl-m  ORAL  PA° T  OF  LOG*. 
“  YSI7f=?.. 

U  n=n‘ 
fPoy=i'  . 

■  xnR!G=YS.I7E  +  i.  ■ 

CALL  ORIGIN' IXORIG.O.r 


A. 

W: 


f  A’Lt  l  ETTF<M0'.?6,0.,0.15,'REAL  PART  OF  LOG- %  90.  »  17* 


OF  LOG. *,90.,  1 
«  YMaX,  Y.M  IN*  IN* 


-'CALI  PLOT  A  (WOPK,  I £Np«  DELX,  YS I ZE,  XLEN,  YMAX.YM IN,  IN,  I B.OX,  INC ,9T)  .  ^ 
P'MAX=Y»»AX  •'  .  **  - 

'pmin=YM!N  "  ’  " 

.WRITEfft.lM  TNDEXP.RMAX  ^  '  .  t 

14  FORM  AT  {//.IX,*  INDEX  OF  PEAK  =• U 10, LOX, ‘MAX.  VALUE,  REAL  PART’ OF  L 
10G.  1PE19-7)  >  ;  .  -  -l  ’  *•  /  v  • 

‘FIND  .INDEXES  •ABi^tft^JHE  PEAK  ’FOR  WHICH. VALUES  RF*W£EN  ARE  ABOVE 
oTHF  CL  IP  IFVEL,  '  .  .  .  .it.' 

INOFXn=I  END  ,  .  -■  .  .  ,  •  / 

INDFXA= 1  .  '  ..  . 

1 F ( OBLT M.L T-0. J  GO  TO  ZO  '*  ,  - 

•CL  I  P=f  l  -  /  ( 10**  ( DBL  I M/,20. 1 )  I  *  .  ‘  -  ‘  * 

.  CL  I  P=XLOG( Ct  IP )  ‘ 

ACL  IPrABSICLTPl  i'  « 

Cl  IP=RMAX*CL TP  '  '  .  " 

IF(  INOEXP’.EOil)  GO  TO  IT  .  „  ' 

NTIHFS  =  I NDF XP- 1  *  ! 

'  DO  14  J  =  1,NTJMES  .  , 

0  K=  I NDF  XP-J  ■  •  .  •  •  .. 

■  indfxa=k  *  .  ’• 

IF{  f  RM'AX-WORK(K) ) -GE.ACL IP)  (jO  TO  .17  , 

:1 6  CONTINUE  ‘ 

’1  7  CONTINUE  «  .  ’  .  •  J  . 

IFCIND6XP-.  EQ.IFND)  GO  TO  20  '  ■  „ 

MTIMFS-INDFXPf 1  1  -  •  '  • 

*  OD  1 9  J=NT I MgS , I ENn  •  '  4 

.  „  I F  (  I  pm  a  X-WOR  K I J  1 GE- ACL  I  PL  ’GO  TO  20  % 

INDFXB  =  J..  .  ^  . 

1 9  CONTINUE  .  ‘  '  » 

70  CONTTNUF  \  .  ,  , 

DP-  77  J=  1-,  I  END  /  ■  ■  -  .  ;*  "  * 

7?  WDR<K(J1=AIMAG(X(J)  1  '•  .  - 

GO-  TO  (610,420) , TA ARON'  ’ 

610  CONTINUE  ° 

CALL  AARON(X, WORK, COR, INDEXAtlNDEXB, CLIP.NPTS, lENO.PI.TWOPI) 

.  GO  TO  630 

470  CONTINUE  « 

CALI SCHAFP (X, WORK, COR, I NOEXA, I NOEXB.CLI P, NPTS, I  END, PI ,TWOP I » TFMP ) 
6^0  C.DNT  TNUF  „ 

WR I T  E ( 6, 7 1  )  INOFXA, INDEXB 

71  FORMAT  (//.  IX,  •  INOEXA  =  *  ,  1 10 , 1  OX,  MNDE  XB  =*,110). 

ROTATE  THE  PHASE  .  " 

TFCINOFXB.FO. IFNO)  I NDFXB= I NDEXB- 1  -  '  •  a 

*  'J^FI  T=  I NDFXB-.I NDFX  A ♦! 


services  group 


'  .\  • 


00  600  J  =  1  .NPFTT 

600  TEMPI v-  ^ 

°SLnPF=WORK(INpEXR)/TE»»P(NPFITl,  1 

K=l.  b  *  o  ’  ‘  -  • 

HQ- 601  J=r»10FXA,INnPXB  '.  /  -  ,  * 

SAVFC J J  =  SLOPF*T-F«p(K)  *  ‘ 

TFMPJ  J  1=W0RK l J 1— S AVEf  J) 

601  K=K*1  '  *  •,  \  . 

HO  60?  J=INDEXA, INDEX*  }  '-  '  * 

f  REALP=R£A1  (X{J»|  .  •  ~  ‘ 

.  ■  .  CORf  J»=CGRCJI*SAVECJ)  . '  : 

60?  X(.!»=CMPLXn*.EALP.TFMP(  J>)  ' 

I P?0  f  •  \  •  ■  ■  • '  * . *  .  . 

tCC=NP  TS/J  ■  -  .  ' 

;..00  51  ^?.rrx  .  .  ...  - 

RFALP^REALlxijI  1  ' 

AI=-AIHAG(X{ Jll 

corinpts-ip)=-cor( ji  •*  '  ■' 

.  X(NPTS-fP)=CMPLXCREALP»An 

5 :b.  IP=1P*1  ..'o' 

TF( IAARON.Fd.21  GO  TO  835  -  ’..*••• 

PLOT  REAL  PART  OF  LOG’.  £ILPPE0 
YS1ZE=2.  / 

IN=1  ^  .  .  •  '  :  . 

iRnx=r  . 

YMAX=PMAX 

YMTN=PMTN  *  ' 

XOR I G*YS 1 7  F+ 1 .  *  *  ...  s'  ~? 

CALL  ORIGIN  C XOR  I G * 0  ’  ‘  »  .  (  ‘ 

CALL  LFTTERCO. 25.0., 0.15., 'REAL  PART  OF  LOG.  CLIPPED*  *9Q.  *’?5)  . 
no  54  *j=l  *NPTS  ■  *  ;•  ( 

54"  wdRKfJ)=REALCXUn  . 

CALL  PLOTAf  WORUrNP-TS*  DFLX#  YSIZErXLEN*  YNAX»  YMINt  IN* I  BOX*  INC*DT) 

835  CONTINUE 

PLOT  i mag.  PART  OF  LOG.  CLIPPED  AND  CORRECTED 
YSt7E=2.  ,  v 

TN=0  *  ” 

IR0X=1  .  ‘  •  ■  «*  .  .  '  . 

,  X0RIG=YSIZE+1  .  -  ,  “  'i 

CALI  QRIGINIXORiG.O.V 

CAL!  I  ETTFR(0.25VQ.,0.15;*IMAG.  PART  OF^OG.  CUPPED  ANO  CORREbTED 
1  *  *00. *40  1  'v 

.on  55  j=i,npts 

55  vioRK(Jl=AI«AGIXf Jl.l  - 

CALL  \PLOTA  ( WORK.NPTStOELX.  YS'l  ZE*  XLEN*YWAX*YHIN*  IN*  !BOX  » INC*  DT ) 

SAVE  X  ARRAY  IN  SAVE  <  '  ' 

00  200  J=r.NPTS  • 

700  SAVE! J )  =  XC  J 1  ; 

TAKE  1 FFT  OF  ARRAY  X  .  * 

CALL  NLOGNIX.NPTS.-^l.OI  ^ 

PI.OT  REAL  PART  *0F  IFFT  (COMPLEX  CEPSTRUM1 
YST  ZE=2.  *  * 

IN»0 

I BPX=1  ...  g, 


,  services  group 


y-srv .  •. 


i  c  " 

PART 


OF  IFFT  (COMPLEX'  CEPSTRIJ.M)  •  , 


r  ^ 


*n 


3i 


201 


Xno  T  C= Y  5  T  7  F*  I  •  ' 

CALL'  OP  1 01 M  (  XPR  IG»  0.  ) 

CALL'  L  5TTFP  f  0i?5,*O»., 0.  15  ,  *RFAL 
1 90.  ,35 )  ,  *•  •  ^  * 

PLOT  C£P%  TP1J"  FRdM  '-W+l  ■  TO  +W 

jdts^nrts*i  ,  ■ 

NHAIf=MPT  R/?‘ 

on  ’57  ‘.1  =  1,  tENO  '  4*  - 

M0P<(J*NHALF-1  1-RFAL.(XfJ*n  •  •  _  *.  * 

■5T OP  ( J  +NHALF-1')=W0.RK  LJ+NWAL  F— 1)  J  -  ‘ 

NHALF=0MALF-1  .  V  *'cr  •  . 

00  -31.  J=!  .VhALF  ,  *  *  '  %  -  .s*  ■  ,  ■ 

wnpkf  j)  =pr ai  (or ENO+j} )  ,v  ■  • 

,STORCJ)=wdRK(J)  S'  .  .  ,  '  *  •  •  •  ;  . 

CAL!  pf  OTA  (  WORK  Vji^TS  »DFL  X,  YSi  ZE  »  X*LEN«  YMAX *YMI!il*  I N*  IROX *  INC*DT ) 
TALI  WRCP(  IOOUT  •fviri.RK  ,NPTS»I » 1 » ID! S»  IDA TE* NAME »  AUMJ 
inonT- 1  noiir 1  "  '  '  " 

on  .201  j=i.,Kipt5  * 

PFA!  °'=PFAL{5AVE(  J)'l  ‘  - 

Xf  ,JI  =ri“Pl  .X!REALP..O,.’j  "  /  ’ 

.  SU«P=YMAX*10.  *  . 

•  SllMN=YMfN/10.  •  ^  '  '  .  -  •  >  ' 

no  49?  f»l,?  '  ,  ■  •  •  '  .  .  . 

DP  493  j=l,NPT5  ° 

!Ff<;rpp(jn  494,403.495 

IFf  STOP (  J I .L T. 5UHNI  STORIJI^SUMN  1  ' 

GO.  TO  493  ..  '  •  . 

IF(STORCJ)  ,GT. SUMP }  <4iT0R  ( J l=SUMP  ”  .  ' 

CONTINUE  *  '  . 

0vST  7E-? .  .  .*  /■  '  ' 

[  N=n  /  ..  . .  ,  •.  • ' v  ' 

IRPX=1  '  .  .  •  . 

X0RTG=V5IZF+1 J 

"CALL  ORLGTNfXORICjO.J  .  • 

PALI  PL  OTA  l  STOP  ,NPT^,06LX,YSirEt  XLEN,  YMAX  ,YMIN*  $  N,  I  BOX*  INC.DT) 

SU'FPsSlJMP/IO.'  '  *  *  . 

su^N=<;u«M/.i6.  '*  .  x. 

continue  - 

i  CALL-  NL-OGN  IX.NpT5,»1.0)  ' 

PLOT  REAL  PART  OF  IFFT  OF  REAL  PART 

YSI7F=?.  „  .  ,  i  *.  •  S 

'■  TN=0 

TPOX=l  ,  f 

X0RIG=VSI7E+! .  ’  '  , 

'  CALL  nRIGTKKXORIG.O.I 

«"Al  I  L  ETTcr>  (  ^.PR.O.  ,0.  l5,*R%P«tFFT.R  .  °*-'»90»  ,  L2 ),  .  ‘ 

00  203  J=!iNPTS  ,  *  .  *  .  *’• 

TFM^J)  =  PF,At  (X(JH  *  .  "  •  /  , 

ST0R(J)rTc9Plj)  °  •  ' 

203  WORK (J)=RFAL (X(J) )  -  .  ‘  - 

CALL  PLOT  A  (WORK  .  NPf  S  »  QEL  X,  VS  I  ZE»  XLFN.YMAX,  YMIN,  IN,  IBOX,  INC^DTI'.. 
FXPANDFO  PLOT  OF  R,\P.  f  FFT.R-.P. 

<;UMp=YMAX/loi  . 

S0mV=Y»>!*N/10.'c3  '  .  ‘  - 


■  t 


494 

495 
493 


40? 


\  r« 


L 


■;  L 


:-1- 

t 


; )  4 


L  ; 
L  ? 


L 

t 
'  i 


T  ” 

’  .  I 


A-  10 


services  group 


l 


nr  496  1=1,?  # 

no  004  j.=  l,NPT< 

IFfSTORCJll  903,904,905 
90^  fF(STOO  (Jl.LT.SLPINJ  STOR(J)=SUMN 
.  GO  TO  904  ,3  '  *  -  * 

905  IF(STnR{J#^f  STOR { J 1=SUMP 

-90.4  CONTINUE.  . 

YST7F=2.  , 

TW=0  •  -  . 

V  I  BOX=l.  u  •’ 

XOR  ro=.Y5T7F+l. 

T  Alt'  .OR  f  GI.N CvXCR I Gi  0.  ) 


o  B 


i  on  I  *  &  .  »*  r  \ 

CALL  PL  CTA'(STOft.,NPTW)'ELX,YSIZE,XLEN,  Y«AX*Y#TN,.lN,  I  BOX,  TNC,  DT) 
•<U5fP=SUMf>/10,.«  *.  .  '  -V  'w  •  X  •*  »,  *• 

5i;mn=shmn/io,  '  '*  ,  °  -  ■  .  .  >  ■  ‘ 


t 


c 


( 


406  CnNToIWUE  f  ■  ■  '■  '  <  - 

‘  DO  205  J  =  r|NPTS  ' 

AI  =  ATMAG(>SAVE(  Jl  >\  .  . 

205  X.(J|=r»PLXCO*.,/'HI 

CALI  NLOGN  (X,NPTS,*i.b»  %  '  '\;l 
*  PLOT  REAL  PART  OF  IFFT  OF,  f MAG  PART 
Ysize=?:  *  -  •  s •  v 

^  <TN=0  .  - 

IBOX=T  .'/o'  •  .1  '  .  . 

XOR  TG= YS I 7E+L  .  0  .  *\  «v  •  *  .  - 

^CALL  OR  IGTNI  XOR  1G,  0.  )•  ■  j'  <, '  ••  <?• 

CALL  LETTER (0.25, 0.70.15, «R.P. I FFTi I.?*,, 4jjj,  1*1 
00  207  J=l,NPTS  '  .■  • 

WORK(J)=REALt(X( J)  J  ' 

207.  STOP  (J)=WORK;U)  *  ..  .  *  «  ‘  ' 

.  fALL '  PL  OTA  (WORK,  NPT  S,  OFLX,  YS  IZE,Xl  EN,  YMAX,  Y.M  IN^ IN, XBOX,  INC»DT) 
EXPANDED  P?.OT  OP  R'.P.  IFFT.  I.P.  "  '  i  .  •/ 

S{»MP,=YMAX/i,0.(<  •  .  • 

’  SUMN=YMIN/10.°  ‘  „  <  .*  \  « 

DO  497  1=1,2 
'  .  4)0  914’  J=l,NPf$ 

IFtSTORCJn  913, 

913v  1F(ST0R(JI.LT.SUMN)  S.TOR  (  J  )=SUMN 
GO  TO  914  ■'  . 

915  IFCSTOfif  J).;GT.5UMP1  STOR ( J)=SUMp, 

914  CONTINUE  '*  *  r  ,  V — /  •  *•  •  ,  v  •  . 

'  ,  .  YST  7E=2  •  •  .  *  \  '  ‘  •  .  v  ^ 

IN=0  *  j,'  *  »  - 

0  ,TR0X=1  ‘  ‘  „  ’  , »,  :  9"  ‘  ' 

XORlG^YStZF+l.  A'”/  •«  ^  V 

.CALL  OR  JGTNC'XOR  IG,  0. 1  f  '*  ■ 

CAL  L  "PLOTA  (  S  TOP,,  NPT5,  DELAYS  1 7.  E  ,’XL EN.,  YMAX,  YMIN-«I  NiVftOX  ,  t  NC  «  D.T ) 
SUMP*SUMP/.ld.  *  .  '  "  ‘  ,  •  -6  ^  ,  '~a 

,sumn=summ/io. 

497  CONTINUE  ‘  '  .-4  *  * 

CALL  ®TT ASK (T-2 1  *  •  \  *  * 

'  WRT  TF(  6,88881  .T7-  ‘  «  *  •  ” 

RR8B  PQRYATC//, IX, ‘TIME -PART  ONE  . * • ,F10.0, »  SECONDS’  1 
. IF( IFILT.F0.01  GO  TO  999  ’  •  •  ‘ 


'•  u 

v  '  ^  •  O 

►  *  *'  :  .  x  »  .  *  t  y  «  •  • 

1,914,915.  *•;  ;  1  \  l 


r 


i 


■j 


-A- IP 


services  group 


<  r 


CALI  <-.TASK' 

nr  Aa  .i=ii^or<; 

SAVFf J  *=rMOL1f  I  TF**»  f  J1  «UfQ9K(  J)  ) 

RFAO  riLTFR  OATA*  - 

rn  r«  t^st  rco^pim, 
nr  ?^rILT=l*  IFIlT 
on  3-oo 

TF*«PfJ)=PGALf<;AVEf.Jll 
MnstKc.n=Ai«AGi«;Avefjn  .  * 
opAm^.m?!  !TYPF*iniTS»T,mx,JTAPER 
rnpv&u'tliQ) 

W'P  ITFC<.-»?n»  .A 

FORMAT  I  //,  l  v,  »  I T Yoc  ,  mi  tpt ,  I D* ,  JT&FR *  1 
SW'lTpf  6.3*1  I  1TYPF.I*ITPT.U2£.JTAPER. 

fopmati  lx.Aum  *-  * 

IFIN»*T=i*iOTS-f?IITPT*2 

CO  TO  f 306.310.3151. ITYPE 

shopt  pas<;  sy<;tf« 

•OP  306  .l=r*4lTPT,  IFIKPT  * 

WPPKf J)=0, 

TF«P!J!=0. 

IFf  JTAPEO-FO^l  I  CO  TP  31'6 
lTAPFP=fNITPT-4  ? 

orc=o.  '  .  .  ! 

no  A  AS  .1=1  TAPFP.INI  TPT 
AQG=DFC/57.79S8 
ARGsfl.4-C.0Sf  ARCH/?. 

WOPK  f  J  )=MOP.K  (  J  1*A0G 
TF^Pf  JWTFYPf  J)*ARG  ; 
nFC^nrG+65.  ,  - 

I T  \  P*f°  =  I F IKP  T  *4 

nec=i80.  •  -t- 

DO  667  J=I F  T  N°T»  I  TAPER** 

ARG=DEG/57.?958 
APC= (1 « •*-CPS(ADGll/2« 

WPPK(  JI*4PG 

TFMP(.n  =  TF''PfJ)*A«G  >*  «. 

.OFG-nFG+45  a  ) 

CO  TO  316  .  | 

LONG  °ASS  SYSTF« 

OO  311  J=! « I NITPT  ‘  '  .  '  ' 

wn<?K(J)=n.  -  - 

tf«pfvn.=n. 

no  312  J= I F I N°T « NpTS  *  J 

WORK  f  J  )  =0. 

T^F^P-f  J  1=0.  '  _  -  ' 

1 F { J  T  A  P  F P  . c0 .1)  GO  TO  316 

ITAPFP=IMITPT+4 

Of  G= 180. 

no  666  •/J=I,NlTPf  ,  ITAPFR 
ARG=H5 d/57, 2958  .  ' 

ARG=(1  l.+  CrS(A»GM/2.  . 

WOPK(J)=WnRl<{J)<=ARG  J  ■* 

.Tf’PP  (.11  =^TF  MP  <  J  )  *  ARG  J 


-\  •  ' 


666  DEG=DFG+45«  / 

TT6PFR=IFI^PTi-‘ 

’  DEG=0.  *'  *  '  *  -  .  '■ 

-t  ’  DO  668  J=fTAPFR.IFlNPT  I. 

A»G=OFG/57.2958'  -  '  • r. 

,  ARG=fl-*COSJAonn/2. 

ar  HORKCJ)sWr.aiC(jJ*A«G.  '  •  ' 

TF*P{JJ=TFPP{Jf *ARG‘ 

668  DEG=DEG*45.  *  * 

s  GO  TO  316 

r  \T.n»«8  riLTFR  ( I KI T PT=NUMBJER  OF  NOTCHES J 
316CONTTNUF  •'  ' 

RFAOC5.31SJ  M NOTCH ( J  J  f J=1 » I NITPT ) ' 

318  FORHATC 81101  *,  " 

Do(-3?2  J=t*INITf>T  “  -  ’ 


l 


£ 


inFri=  inotchc  Ji  -  mx 

.JOF#?=  TNOTCH (  J  l  +  inx 


K 


I ! 


c 


Pa 


f; 


E 


DO  323  K=IOFX  J*  IOEXrz  * 

WORK (K 1 -0» 

3?3  'TFMPCK  1  =0»  .  v 
322  CONTTNUF 

.  IF.f  JTAP.ER.EO.il  GO  .TO  316 

316  CONTINUE  ~  . 

DO  317  J=l,NPTS 

;  .  *  REAL°=TFM»C  J  J  ^WTJRKC  J)  ^ 

317  XC  J)=CMPl  X ( RFAlPiO^l 

take’  fft  of  array  X 
CAL’i.  NL0GN<X.NPTS,-1.0)  * 

REVERSF  ROTATION  ~  , 

DO  *»9  J=1  •  NPTS 
RFAl.P-PEAi  (X(  J.J.)  ■ 

•  Af=AlMAG(X{  J.H+CORI  J) 

69  X(  J1.=C^PLX(RF*LP,AI1 

Pl_OT  RFA^'PART  OF  SECOND  FFT 
,  ,_.YS.TZE=2«  •  *  '  -? 

IB0X=1  -  • 

YHAX=PMAX 
YMTN=PMIN  . 

X0RIG=YST7F+1. 

/  p AL 1  OR!GIN(XnRIG»a.K  •  •  ,  .  ‘  *  *  * 

,  CALL  ,LETTER(0'.25.0.f  0.1.5* ’REAL  PART  OF  SECOND  FFT  *  ,90.',  23  ) 
on  64  j=i,.npt.s  '  '  -  ,  V  •  . 

64  HORKCJ)  =RFAL(XC.JH  .  ’  "  :  '  '  .  ,  .  - 

?CALL  Pl.ptA  (^ORX*  NPTS* DEL X*  YST  ZE  »  XLEN*  YHAX  *YMIN»  I  N*T  BOX  »  TNC"V  OT 1 

#  TAKE  EXPONENTIAL"  OF  ARRAY  X  ,  *  \  • 

<9  \  y 


.  / 
**/ 


✓ 


f? 

a 


Ml. 


DO  60  J=1,'NP^S 
60  X'(  J  |'=C EXP IX1JII 
-  .  *  I F (  TTYPE .EO .2  )  ;G0  TO  450 

GO  TO  <  I50,'l60)!,  rsw 
X  160  DO  170  J-KNPfS 

REALP=-REAL1X(J)  1  •" 

AI=.-A'T  MAGFXI  J)  )  • 

•  170  X(  J1=CMPLX(RF6LP»AI) 


v. 


v 
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ISO  CO'JTINUF  - 

/F  ’  ,  TAKF  IFFT  OF  AOS  AY  X  '  ' 

:  CAM.  N(  OGN(X,NPTS.«-1.0)  .  . 

\c  PI  OT  RFAL  PART  OF  SECONO  If  FT  (OUTPUT  I 

.  "Y5I7F-2.  .  *  -'•  '  *  .  .  .  •  * 

*  -im=  i  *  •  .  * 

I  a0X  =  7  -  ,  .  '  ,  -  *  -  ‘  .  ’ 

_  YMAX=f  MAX'"- — -  '  -  .  •  \  '* 

•YMTN=TNIN  .  ->  ... 

XORIG^YSIZE+l .  •*’  *'  *  -  .  ; 

CALL  ORIGIN  ( XORIG .0. )  *  ,  . 

GO  TO  (791  *752 •  79,3)-,  1  TYPE. 

j  *  701  CALL.  LFTT/FR(0.?5,0.,b.'15i*REAL  PART  OF  SECOND  IFFT  (SHORT  PASS  OUT 
.  .  ~!PUT  )  90.4  44)  •  '  -  -  -  •  '  - 

!  GO-  TO  794-.  ;  '•  ‘  '  _  : 

i  792  CALt  LETTFR(0.2S,0.,0.15.*REAt  PART  OF  SECOND  IFFT  (LONG  PASS  OUT'P 
.  ‘  '  1UT) * .90..43)  *  ‘  ^  y  .  ' 

J*.  GO  TO  794  '  X  ‘  ; 

-  793  CALL  I  ETTF«^0.^5,0.;0.V5.»REAL*  PART  OF  SECOND  IFFT  (COMB  OUTPUT)*, 

*•!'.'  190. .38);  \'  -  4  ,-v  /  j 


!  794  CONTTMIJF  *  *  •*■...  . 

DO  62‘.J=1,NPTS  -  -  '.y  .  1  1 

'?  wn»K(  j)=PFALn$m'>  .  n  -•*.  . 

IF(  ITYPE.EQ.2r  GO  TO  $85  \.  '  - 

•DO’  760-  J=  l, 'LPT S  •  ■  .  1  _  '  > 

AJ=— J  ••• 

760' WORK  (J)  =  HORk'(J)>(  ALPHA**AJI  '  - 

GO  TO  6.86  .*  *  ~  ’ 

,>,86  IN=0  .  -  *  '  *  ' 

K^l  -  *  v 

-  DO  687  J=1 ,LPTS  f  *  *  '  .  ,  . 

*  IF  (-WORK!  JI  .GT.D.5)  GO  TO  688  •  ; 

!  6*7  K*K+l  ' 

'  \  K*r  ’  •  ^ 

.6^8  CONTINUE  .  c'  '  •  *  ' 

*  ■  HR  I TFC  6,690)  K  •  ‘  *',.„■ 

,  -690  FORMAT  (7  7. IX^' INITIAL  POINT;  FOR  UNHE IGHTING  LONG  PASS  OUTPUT  =•  , 

\  -  '  II 10)  .  .  -  '  -  ‘ 

l»lpts4k-i  •  ■  /♦.y  '  *"• 

m=o  .  •«  .  - '  /  '  ■ 

00**639  ^J  =  K * L  V  .  r  ' 

;• '  s  •  /  r  '  -  -  •  *•*  ’ 

A,J=-M  .  .  *  ’ 

689  W0RK(  J)  =  HORK(  J)*  (  AL<PHA**AJT  ..... 

6R*6  CONT INIJF  *  .  -  /  ' 

•CALL  PLOTAl  WORK, KPTS,AE  LX.  Y  SIZE*  XLEN»$8AX*YMIN»  IN  ♦’I  BOX  »'INQ»DT )  ■ 

,  .  CALL  WRCD(  IOOUTv  .WORK,  NPTS,  1-,  lflDTSTlDAtE, NAME,  AUH)  • 

-  inmfTr  tooutm  .  /j  ,  \ 

'•  IF(  TTYPF.FQ.2)/gO  -TO  888  1  *#  .  \  . 

f*  PLOT  -Til AC F  (INPUT  SIGNAL  MINUS  OUTPUT)/t  ,  ;\  ^ 

-  TN=1.  '  «  *  /  /  -  ’  *  • 

I90X=>  /  /  */'**•. 

Y  m  A  X  =  T  WA  X  /  '  *  /  '  *  *  , 
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/ 


$ . . 


YMIN=T*TN  j  . 

X0Rt'G«VST7F+l.  ,  1  - 

CALL  ORIGIN  !XORIG*Q,.) 

•  CALL  LFTTERIO.25,0.  ,Ow 15» • INPUT  SIGNAL  MINUS  OUTPUT*  ,90.  ,25 1  : 

00  31  J=1,NPTS  ,*  ‘  c  •  -  ' 

81  TEMP  l  J  )=TR  ACE  f  J  )— WORlfr(-J  )  .  ' 

CALL  P.LOT,A(TEHPtKPTS*  AELX»VSIZE» JCLeN#  YMAX»YMIN.»IN*IBOX,  INC»DT) 

-  CALL-  WRCPC  IOqUT,TEMP,NPTS.l,l*mS,IOATE*NAMEviAUMl 
■  -1  •  !DOUT=IOOUTM  -  <  .  \ 

888  CONTINUE  ’  •  '  . '  .  1 

CALL  TTASKIT2) 

•  WRITE!  6,81199)  IFILT.T2  *  -  -  - 

3889-  -FORMAT'!  /  /  ,1X**TIMF*,I5**  FILTERS  NF10.0,*  SECONDS* ) 

999  COWINUF  -  •  -  ' 

'  CALI  ORIGIN!*.,!).)  -  - 

,  C/LL  BLOKNO  .  .  '  '  l  •  .  . 

1000  OONTTNUE  :  ‘ 

1009  /iONTINUF  '  *  ,  .  . 

.  /call  fndplt  *  • 

/  CALL  NN100*  ' 

\stop  *  ® 

’END  '  (  .  ' 

SUBROUTINE  PLOTA!DATA,*NPTS,DELX,YSIZE,XLEN,YMAX,YMIN,  IN,  I  BOX^  TNC, 

'  IDT)  , 

DIMENSION  DATAU)  *  *  •  •  /  # 

CALL  BLOKNO  '  ' 

XLEN=! NPTS-l )  *OELX*INC*DT  . 

IFIIN.EO.l)  GO  TO  10 
'  YMAX-DAT  All)  . 

YMTN=DAT  All)  . 

DO  1  J=?.»  N°TS  /  . 

IFIDATAljy.GT.YMAX)  YMAX*DATA { J) - 
1  IF1DATA! Jl.LT-YMIN)  YMIN*DATA!J) 

10  SF*YSIZE/!YMAXtYMIN?  '  ♦'s 

IF!  IROX.FO.l  j  CALL  SQUAREI-YSI  ZE»(y.  ,0.  ,XLEN) 

N=6  j  '  '  •  .  •  .1 

no  4/.J=l,7  ",  *  ■  ft 

IF! ARS!YMIN*10**N).LT. 8400000.)  GO  TO  5  V 

•  4  N*Nf  1  ‘  ■  sf  \ 

GO  |T0  6  . ‘  '  \  ,  ' 

5  CALjl.  RIGHJ-J!0.iO.|,Oil,.YMIN,90i,N)"  '  ' 

~  6  I  FljYMrN.GT .0.  .OP .YMAX.LT .0.  )  GO  TO  2  .  A  1  , 

ZERO=YMf  N*SF  .  _ 

S  CAljl  RIGHTJ!ZFRO,0.,0VTTOtO^9Q.jX)__ 

CALL  PLOT! ZERO, 0. , 3 )  "  ", 

'  CALL  D  A  SH!-ZERO*  XLPN ,0.1  ,  G.  i )  „  \  •  C  . 

7  N=tT>  '  -  b  ■ 


DO'  7  J  =  1 , 7 
/{  F  (  ARS  ( YMAX*I  ( 


/lF(A8S(YMAX*10**N).LT  .8400000. )  GO  TO  8 
7/N*N-l' 

/  GO  TO  9 

'*•  CALL  R  IGHTJ(-YSIZE,0.,0'.l,YMAX,90.,N) 

9  SFA=-SF  ' 

AY,=  0. 
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AjUlOA  TAfl)-YMIN)*SEA 
CALI  PLOTIAX.AV.3I 
on  3  J»2,NPTS  ...  * 

AY* t J- 1 1 *DELX*INC*DT 
AX* I  DATA! J)-YMIN1*SFA 
3  CALL  PLOf CAX. AY.2I  . 

•  RETURN  .  .  '  .  *'  -  - 

ENO  '  '  c  ' 

SUBROUTINE  NLOGN  (X, LX, SIGN) 

SUBROUTINE  NLOGN  COMPUTES  THE  DISCRETE  FOURIER  TRANSFORM  -BYS THE 
N  LOG  N  ME  THQO  OR  THE  F^ST  FOUR I£R  TRANSFORM  METHOD 


T*P\ARGUMENTS  ARE  1 

N  •  WHFRE  2**N  IS  THE  NUMBER  OF  .TERMS  IN  THE  X*  ARRAY 
X  THF  ARRAY  OF  COMPLEX  NUMBERS  FOR  BOTH  INPUT  AND  OUTOUT 
SIGN  EITHER  -1.0  OR-  +1.0  .  -  ’  ' 

1  X  COMPUTED  AS  2**N 

NMAX=IARGEST  VALUE  OF  N  TO  BE  PROCESSED 
NONOUMMY  DIMENSION  MfNMAXY  . 

FOR* EXAMPLE. IF  NMAX*25  THEN 

DIMENSION  MI25I  ;  # 

DIMENSION  XIII  }  ■ 

COMPLEJf  X.WK.HOLD.Q'  *  •  ' 

FLX=LX 

ILK*1  ’  '  * 

00  3380  1*1.25 ,  - 

IF  CLX.LE.ILKf  GO  TO  3381  ; 

.  -  N=I 

UK=!LK*42 
CONTINUE 
00  !  1*1, N 
Mm«2**IN-I» 
no  io  L*i.N  • 

NBLOCK*  2** ( L- 1 1 
LBLOCK=LX/NBLOCK 
LBHALF=LBL0CK/2 
K*0 

DO  11  IBLOCK-l.NBLOCR  ” 

FK*K 

V = S !  G  N  *  6 . 2  Q3 1 8 5  3  *  F  K  /  F  L  X 
WK*CMPLX (COS(V).SIN(V) ) 

I START=LBLOCK* ( I  BLOCK- 1 ) 

DO  2  1*1 .LBHALF 


3380 

3381 

1 


/. 


J*I  START4-1 
JH*  JHBHALF 
Q=Xt JH)*WK 
XI  JH»*X(  J)-0 
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services  gro.up 


4 

II 

10  ' 


K=K^H(fn 
CONTTNUF 
CONT  INUF 
K=0 

•  on  7  J= 1 . LX 
IF  (  K  .I  T  .  J  )  GO  TO  5 
H0LJ3=X  ( J T. 

'  X(JI>((KO) 

X(K*11=HOLO 
5  DO  6  l=l,N 

•  1 1  =T 

IF!  KiLT-MUl ) 
ft  K=K-M(  I  f 

7  KiK+Mdt)  ' 

IF f SIGN. LT. 0.0)  RETURN 
no  8' 1=1. LX 

8  XU.I=X(n/FLX 
RFTURN  ° 

FN0  - 

SUBROUTINE 
1  TW0PI1 


GO  TO  7 


AARON (X. WORK, COR. INDEX A, INDEXB.CL IP. NPTS, TEND, PI. 


40 

4? 


GO 

GO 


TO 

TO 


30 

30' 


4,10 


0  KMF  NS  I  ON  X  (1-1  .  MARK  ( 1 1  .  COR  ( l  ) 
COMPLEX  X  . 

I SUB=NPTS/7 

1F( INDEXA.EO.il  GO  TO  40 
WORK (11=0.0 

X ( 1 1=CNPLX (CL  I P  » WORKf 111  % 

MORK (2 1-0.0 

X (2 1 =CMPLX (CLIP .MORK ( 2  1  K 
Gh  .Tf!  42' 

CONTINUE 
CON.TINUF 
ertR(ii=o. 

COR (21=0. 

DO  25  LL  =  3* 1  SUB  ,  . 

IFILL.LE.  I.NDEXA) 

IF  ( LL  «GT .  I.NDEX8  1  ■ 

I F( LL .LF/. INOEXA+2  1  GO  TO. 410 
YHAT=’2*W0RK  (LL-ll-WORK  (LL-2 1 
XMOD=A  MOD.(  YHAT  .TMOPIJ 
I F  (  X  M0  D • GT . R  I  1  •XHOO=XMOO-THOP.l 
IFLXMOD.LE.tPI 1  XMOD=XMOOtTMO>l 
PHA=WORK(LL» 

DTFF=ABS(  X  MOD- PH  A  1  * 

SIGN=+ 1. 

TFIP.HA.  -XMOD, 

I;f  (DIFF.LT.PI )  WORK ( LL  1/^YHAT+PHA 
IF ( D  IFF . GT ,P 1 1  WORK! LL )  =  YHAT+PHA 
r.OR(IL)-=PHA  -WnRKdLf 
RX=REAL('X(LL1  1  / 

X  (  LL  J  =  CMP,t  X  (  RX .WORK (LL 1  1 
GO  TO  75 

CQR(U  )=0.  '  ‘  • 

RX=RFAL(  X(LL)  )•  •  < 


I.L  T.O.  1  SIGN*-!. 


-XMOD 

-XMOD-( TWOPI*SIGN) 


“S - 
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*  X(  LL)-r...MRL^WX,WORK(.LLn 

r.n  T.n  ?s 

30  WORKfLLMO.O  ■  •  ' 

COR  ( LL.)  *0.  ' 

XILL  l=CMPLXrCL  IP»WORK(l.LU 
25  CONTINUE  ’■ 

!F{  INOEXB.LE.  ISUB)  XWEND)=CMPLX(CL  IP, WORK!  IEND)  ) 
COR(  IEN0)=0.  . 

RETURN 

END  „  __ 

SUBROUT  INE-SCIhAFR  ( X ,  WORK , COR NOEX  A*I,NDE£B  ,CL  I P*NPT,S* 
'  1  THOpf  ,'TFfJ.P)  ■..* 

DIMENSION  X(1  J  .WORKfl )  ,C0R(1  I  , TEMPI  1 )  L 

complex  x  n- 

ISUB*NPTS/2  v.  % 

EP=A.8  \  ,  .  ' 

EP-1..  -  \ 

EP^TWOP I-Fp  '  , 

INDEXA=1. 

*  I N0FX8*I END  ' 

THIS  ROUTINE  ALLOWS  NO  CLIPPING  ^ 

COR  III *0.  / 

TFM.PIl  |»0.„  .  S  ■  '  . 

00  1*5  J=  2*  I  SUB 
COR  I j)=CORfJ-lI 
0  IFF*WORI<  f J I  -WORK!  J-?  1 1 
IFCOIFF..GT.EP  I  COR(JI*C0R( J-ll-TWOPJ 
4  IFIOIFF.LT. -EP)  COR  (J)*CORC d-U+TWOPI 
IS  TEMP ( J 1=W0RK ( J ) ♦COR ( J I 
00  20  J~ l t I SUB 

.  COR ( J) =-COR { J )  '  ^  •  . 

RFALP*REAL(XIJ)I  '  > 

WORK! J) = TEMPI Jl 

20  X(  J|*CMPLX(PEALP.TEMP(  J*  I  ;  ^ 

COR ( IEND I =0.  •  ' 

RETURN 

ENO  -  ' 


IEND*  PI * 


1 


